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FOREWORD 

Modem  themes  of  static  electric,  magnetic,  and  gravitational  potential-fields  are  based  upon  the 
early  investigations  of  Sir  William  Gilbert  (1600)  and  of  Sir  Isaac  Newton  (1686)  and  upon 
subsequent  investigations  covering  more  than  two  centuries  by  numerous  other  well-known 
mathematicians  and  natural  philosophers,  culminating  in  the  synthesis  of  electricity  and 
magnetism  by  Sir  James  Clerk  Maxwell  (1864)  and  the  synthesis  of  gravitation  and  geometry  by 
Albert  Einstein  (1916).  Applications  of  electromagnetic  field  theory  and  gravitational  field 
theory  to  the  remote  sensing  of  geophysical  parameters  through  inverse  modeling  techniques  has 
been  a  major  theme  throughout  the  remainder  of  the  2()th  century.  The  Navy’s  interest  in  these 
subjects  stems  from  its  Anti-Submarine- Warfare  mission,  from  its  concern  for  locating  undersea 
Hazards-to-Navigation,  and  from  its  interest  in  Attitude-Heading-Reference  systems. 

This  report  reviews  and  synthesizes,  and  in  so  doing  self-consistently  unifies  through  tensor 
analysis,  the  subject  of  static  geopotential-field  theory  with  respect  to  rectangular-prismatic 
bodies,  which  may  be  considered  as  the  basic  building  blocks  of  natural  geophysical  phenomena. 
Although  this  report  originated  from  geophysical  considerations,  the  mathematical  expositions  it 
contains  are  generic  in  the  sense  that  they  are  also  applicable  to  renx>te  sensing  problems  in  the 
fields  of  classical  physics,  biophysics,  and  engineering. 

D.  J.  WHITFORD 
Commander,  U.S.  Navy 
Commanding  Officer 
Acting 
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ABSTRACT 


We  review  from  first  principles  the  forward  and  inverse  geopotential  modeling  problem  with 
respect  to  the  magnetic,  electric,  and  gravitational  fields.  Under  the  assumption  that  any  region 

of  the  Earth's  crust  can  be  conveniently  subdivided  into  rectangular  prisms  composed  of 

\ 

uniformly  distributed  geophysical  parWneters  (i.e.,  magnetization,  polarization,  density,  etc.),  a 

\ 

wide  variety  of  generalized  Poisson  relations  among  the  three  potentials  and  their  corresponding 

\ 

vector  and  gradient  fields  as  well  as  their  co^esponding  geophysical  parameters  is  established. 

\ 

It  is  shown  that  when  both  vector  and  gkdient  tensor  components  of  a  potential  field  are 


simultaneously  available,  inverse  problems,  iiuch  as  determining  the  depth  to  the  oceanic 
magnetic  basement,  can  be  split  into  a  purely  "geometric"  problem,  which  seeks  to  determine  the 
dimension  and  position  parameters  of  one  or  mor^  prisms  and  a  purely  "geophysical"  problem, 
which  seeks  to  determine  the  physical  properties  (lye.,  the  magnetization,  polarization,  density, 
etc.)  of  the  prism.  The  geometrical  problem  is  nonlinear  and  must  be  solved  iteratively  using 
standard  stochastic  or  least-squares  inversion  techniques,  while  the  geophysical  problem  is  linear 
and  may  be  solved  by  direct  inversion  once  the  geom(  trie  parameters  have  been  established.  By 
splitting  the  inverse  problem  into  two  parts  in  this  mt  nner,  the  more  troublesome  aspects  of  the 
well-known  non-uniqueness  problem  associated  witp  geopotential  inversions  are  minimized 
relative  to  other  inverse  modeling  methods  that  attempt  to  solve  for  the  geometric  and 
geophysical  parameters  simultaneously  using  scalar  or  vector  data  alone.  This  approach  to 
inverse  modeling  does  not  require  direct  measurement  of  the  gradient  field  since  any 
2-dimensional  regional  survey  of  just  one  vector  field  component  can,  through  rectangular 
harmonic  analysis,  yield  a  potential  function,  which  in  turn  can  be  used  to  compute  all  vector 
and  gradient  tensor  components  of  the  measured  field. 


V 


1.  INTRODUCTION 


The  Earth's  crust  is  generally  composed  of  a  variety  of  inhomogeneities,  the  geometric 
parameters  (i.e.,  dimensions,  orientation,  distribution,  etc.)  and  the  geophysical  parameters  (i.e., 
magnetization,  polarization,  density,  etc.)  of  which  we  seek  to  determine  from  scalar,  vector,  and 
gradient  tensor  measurements  of  one  or  more  geopotential  fields  (i.e.,  magnetic,  electric,  gravity, 
etc.).  As  is  well-known,  the  geopotential  inversion  problem  is  inherently  non-unique  in  the 
sense  that  for  a  specified  potential  field  measurement,  an  arbitrary  choice  of  the  geometric 
parameters  defining  the  source  of  the  potential  field  has  a  corresponding  set  of  geophysical 
parameters  that  will  leave  the  observed  potential  field  unchanged.  Alternatively,  an  arbitrary 
choice  of  the  geophysical  parameters  has  a  corresponding  set  of  geometric  parameters  that  again 
will  leave  the  observed  potential  field  unchanged.  Fortunately,  many  of  the  infinite  variety  of 
possible  source  geometry  and  geophysical  parameter  combinations  that  could  produce  a 
measured  field  are  physically  unrealistic.  Unfortunately,  many  other  possibilities  are  reasonable. 
Selecting  among  all  the  reasonable  solutions,  that  one  which  is  actually  the  source  of  the 
observed  field  is  rarely,  if  ever,  achieved.  Instead,  using  supplemental  information  and 
optimization  techniques,  the  "most  likely"  solution  is  generally  considered  acceptable. 
Historically,  however,  only  the  scalar  and  vector  component  measurements  of  a  potential  field 
have  been  used  to  determine  the  geometric  and  geophysical  source  parameters.  Gravity  and 
geomagnetic  data  processing  and  interpretation  techniques  used  during  the  decade  of  the  1960's 
are  summarized  by  Grant  (1972).  Inverse  methods  using  the  scalar  field  have  been  put  forth  by 
Bott  (1959,  1963,  and  1967),  Talwani  (1965),  Richards  et  al.  (1967),  Bhattacharyya  and  Navolio 
(1975),  Bhattacharyya  and  Chan  (1977),  Pedersen  (1977),  and  Bhattacharyya  (1980),  while 
those  using  vector  fields  have  been  based  on  techniques  developed  by  Plouff  (1976)  and  Barnett 


(1976),  among  others.  These  methods  generally  require  a  simultaneous  solutions  for  both 
the  geophysical  and  geometric  parameters  of  the  inverse  problem,  a  process  which 
leads  to  certain  well-known  ambiguities  as  discussed  for  instance  by  Skeels  (1947)  and  Roy 
(1962). 

The  inversion  technique  developed  in  this  report  takes  further  advantage  of  the  gradient 
tensor  field,  which  is  presumed  to  be  either  measured  or  derived  from  a  rectangular  harmonic 
potential.  The  potential  in  turn  is  generated  from  survey  measurements  of  its  corresponding 
vector  field  which  covers  the  2-dimensional  region  of  interest.  The  addition  of  gradient  tensor 
information  permits  the  separation  of  the  geometric  and  geophysical  portions  of  the  problem. 
This  improves  the  stability  of  the  inversion  process.  Advantage  is  further  taken  of  the  simplicity 
of  the  rectangular  prism  geometry  and  of  the  assumption  of  uniformly  distributed  geophysical 
parameters  within  each  prism,  which  permits  the  closed  form  (i.e.,  non-numerical)  evaluation  of 
otherwise  complicated  integrals.  Any  geological  region  can  be  reduced  to  a  collection  of  such 
prisms. 

As  a  prelude  to  the  inverse  problem,  it  is  first  necessary  to  solve  the  forward  geopotential 
modeling  problem  corresponding  to  that  of  a  geophysical  parameter  uniformly  distributed  within 
a  prism.  This  is  an  old  problem  which  dates  back  almost  200  years  and  which  has  been  solved, 
forgotten,  and  resolved  in  piecemeal  fashion  over  that  time  by  a  long  list  of  researchers,  perhaps 
the  first  of  whom  was  Poisson  himself,  who  first  formulated  (1826)  the  connection  between  the 
magnetic  and  gravity  potentials,  now  referred  to  as  Poisson's  relation.  Modem  interest  in 
potential  field  modeling  based  on  simple  geometric  structures  resurged  during  the  early  part  of 
this  century  with  the  works  of  Kellog  ( 1 929)  and  Barton  ( 1 929),  which  laid  the  foundation  for 
the  analytic  results  concerning  rectangular  prisms  for  vector  and/or  gradient  tensor  components 
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of  the  magnetic  and/or  gravity  fields  by  Nettieton  (1947),  Sorokin  (1951),  Haaz  (1953), 
Bhattacharyya  (1964),  Nagy  (1966),  Goodacre  (1973),  Plouff  (1975a,  1975b),  and  Okabe  (1979 
and  1982),  among  many  others.  These  efforts,  however,  did  not  fully  exploit  the  high  degree  of 
unification  that  exists  between  the  gravity  and  electromagnetic  fields.  This,  we  attempt  to  do  in 
a  methodical  fashion.  In  doing  so  and  to  remain  self-contained  and  internally  consistent  from  a 
notational  point  of  view,  many  well-known  results  will  be  rederived,  others  will  be  generalized, 
and  new  results  generated  to  fill  in  a  few  gaps  that  still  remain  in  the  forward  problem  of 
modeling  geopotential  fields  generated  from  rectangular  prisms  uniformly  filled  with  an 
appropriate  geophysical  parameter.  Along  the  way,  the  usefulness  of  rotational  invariants  that 
may  be  constructed  from  the  various  geopotential-field  vector  and  gradient-tensor  components 
will  be  pointed  out,  and  generalized  Poisson  relations  will  be  constructed  and  tabulated. 

Also,  through  the  generalized  Poisson  relations  that  are  derived  from  the  forward  modeling 
results,  the  connection  between  the  potential  fields  and  the  Riemann  curvature  tensor  in  the 
weak  gravity-field  limit  of  General  Relativity  is  made,  along  with  the  observation  that  any 
classical  unified  field  theory  (i.e.,  any  extension  of  General  Relativity)  must  reduce  in  the  weak 
gravitational/electromagnetic  field  limit  to  the  Newtonian/Maxwell  geopotential  fields  of  the 
forward  model.  This  in  turn  necessarily  places  a  constraint  on,  and  provides  a  possible  starting 
point  for,  the  development  of  a  classical  unified  field  theory. 

In  the  discussions  that  follow,  a  local  Cartesian  coordinate  system  will  be  used  such  that  the 
X-axis  is  oriented  toward  the  North,  the  Y-axis  is  oriented  East,  and  the  Z-axis  is  oriented 
vertically  down  into  the  Earth.  The  observed  potential  field  data  sets  derived  from  geophysical 
surveys,  which  are  originally  collected  in  geodetic  coordinates,  must  therefore  be  appropriately 
preprocessed  to  isolate  the  crustal  contribution  to  the  fields  as  they  appear  in  the  local  Cartesian 
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frame.  Our  procedure  is  first  to  remove  temporal  variations,  which  in  the  magnetic  case 
amounts  to  removing  such  effects  as  the  Solar  Quiet  Daily  Variation  (Sq  DV)  and  its 
corresponding  induction  field,  using  nearby  magnetic  observatory  data,  and  subsequent  removal 
of  the  Earth's  core-generated,  long-wavelength  features  from  the  data,  using  an  appropriate 
spherical  harmonic  model  (e.g.,  WMM-90  to  degree  12  for  magnetic  field  data  and  WGS-84  to 
degree  18  for  gravity  data).  Second,  a  rotation  of  the  residual  potential  field  vector  and  gradient 
components  from  geodetic  to  rectangular  coordinates  is  performed  using  the  spherical  coordinate 
system  as  an  intermediate  stage.  The  corresponding  coordinate  transformations  are  also 
performed.  Finally,  biases  and  linear  regional  trends  are  removed  from  the  residual  field 
components.  The  resulting  residuals  of  each  field  component  are  then  assumed  to  be  generated 
by  the  geophysical  properties  of  the  local  crustal  materials.  It  is  with  these  properties  that  we  are 
concerned. 

2.  THE  MAGNETIC  FIELD 
2.1  Theoretical  Background 

Because  the  electromagnetic  field  obeys  the  superposition  principle,  modeling  the  magnetic 
field  due  to  non-uniformly  magnetized  oceanic  crust  in  some  region  of  the  Earth,  which  may  be 
characterized  as  a  collection  of  uniformly  magnetized  prisms,  reduces  to  the  analysis  of  only 
one  prism.  This  analysis  begins  as  usual  with  Maxwell's  equations,  which  in  Gaussian  units  are 
given  by  the  following  relations: 

V*D  =  4n(j  (la) 

V.B  =  0  (lb) 
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VxE  +  x  =  0 

(Ic) 

VxH  = 

(Id) 

where  all  of  the  parameters  have  their  usual  meanings 

constitutive  relations: 

as  do  those  in  the  accompanying 

D  =  E  +  4nP 

(2a) 

H  =  B  +  4nM 

(2b) 

and  the  generalized  form  of  Ohm's  law: 

J  =  -n*E  (3) 

where  i\  is  the  conductivity  tensor.  A  dimensional  analysis  of  the  constitutive  relations  shows 
that  the  units  of  the  electric  displacement  D,  the  electric  field  E,  and  the  electric  polarization  P 
are  equivalent  in  Gaussian  units,  which  are  taken  to  be  statvolts/cm.  Likewise,  the  respective 
units  of  the  magnetic  induction  B,  the  magnetic  field  H,  and  the  magnetization  M  in  Gaussian 
units  are  also  equivalent  and  are  taken  to  be  nanoTeslas  (nT).  The  electric  charge  density  a  has 
units  of  statco^llombs/cm^  the  current  density  J  has  units  of  statamps/cm^ ,  and  the  conductivity 
tensor  i)  has  units  of  sec  '.  The  speed  of  light  c  has  units  of  cm/sec. 

In  the  discussions  that  follow,  a  local  Cartesian  coordinate  system  is  used,  which  is  oriented 
so  that  the  positive  X-axis  is  directed  North,  the  positive  Y-axis  is  directed  East,  and  the  positive 
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Z-axis  is  directed  vertically  down  into  the  Earth.  It  is  assumed  that  the  field  observations  have 
been  reduced  as  explained  in  the  introduction  and  that  we  are  dealing  with  the  resulting 
residuals  that  are  generated  by  the  Earth's  crust  alone.  Then,  under  the  assumption  that  no 
electric  fields,  currents,  or  charges  are  present,  the  magnetic  fields  are  necessarily  time 
independent,  and  Maxwell's  equations  reduce  to  their  magnetostatic  forms: 


V*B  =  0 

(4a) 

VxH  =  0 

(4b) 

with  eq.  (2b)  as  the  only  remaining  constitutive  relation.  The  second  of  these  equations  implies 
that  the  magnetic  field,  H(r),  is  the  gradient  of  some  scalar  magnetic  potential  0(r), 
(units:  nT-cm).  That  is: 


H  =  -Va)  (5) 

By  taking  the  divergence  of  eq.  (2b)  and  combining  the  result  with  eqs.  (4a)  and  (5),  the 
following  familiar  Poisson  equation  is  obtained: 

=  4n  V  •  M  (6) 

It  has  the  well-known  solution: 

W=-f  (7) 

j  V/ 
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Here,  we  use  the  notation  of  Jackson  (1975),  where  the  primed  coordinates  range  over  the 
magnetized  volume  V,  while  the  unprimed  coordinates  refer  to  the  point  of  observation.  This 
expression  for  the  magnetic  potential  is  valid  regardless  of  whether  or  not  the  point  of 
observation  r  is  inside  or  outside  the  magnetized  volume. 

Eq.  (7)  may  be  written  in  more  convenient  form  by  employing  the  following  identity: 

V*(4'M)  =  +  4'V*M  (8) 


where,  if  the  scalar  function  'F  is  identified  as: 


'F(r,r') 


(9) 


then: 


Using  the  divergence  theorem,  it  can  be  shown  that  the  first  integral  vanishes  over  any  volume 
that  completely  encloses  the  volume  V.  In  the  second  integral,  use  may  be  made  of  the  rule: 

Thus,  the  magnetic  potential  reduces  to  its  most  useful  form: 
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(12) 


W  =  -V*  f  ,  dV' 

Now,  measurements  of  the  magnetic  field  are  usually  made  in  regions  of  space  that  are  void 
of  any  magnetic  materials,  as  is  the  case  with  aeromagnetic  surveys,  for  instance.  In  such 
regions  the  magnetization  vector  M(r)  is  zero,  so  that  by  eq.  (2b),  H(r)  =  B(r).  Consequently, 
eq.  (5)  reduces  to: 

B(r)  =  -  V<I>(r)  (13) 

The  gradient  of  B(r)  yields  the  magnetic  gradient  tensor  5B(r)  (units:  nT/cm): 

^(r)  =  VB(r)  (14) 

Both  B(r)  and  ®(r)  are  measurable  quantities.  The  state  of  the  art  in  magnetic  gradiometry  is 

discussed  by  Fram  et  al.  (1974),  Kekelis  (1984),  and  Hastings  et  al.  (1985). 

In  the  local  Cartesian  reference  frame,  the  magnetic  induction  and  its  corresponding  gradient 
tensor,  expressed  in  terms  of  their  respective  vector  and  tensor  components,  are  written  as: 

B(r)  =  B,(r)i  +  By(r)j  +  B.(r)k  (15) 


«(r)  = 


f 

dB, 

dB, 

dx 

ay 

dz 

dBy 

dBy 

dBy 

dx 

dy 

dz 

3Bz 

dBz 

dBz 

dy 

dz 

(16) 


I 

I 

I 
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Rather  than  using  the  above  notation  to  describe  the  vector  and  gradient  tensor  components  of 
the  magnetic  field,  we  will  find  it  more  convenient  to  use  tensor  notation.  In  tensor  notation,  the 
coordinates  and  conesponding  subscripts  x,  y,  and  z  are  replaced  by  numerical  indices  ranging 

from  1  to  3  such  that  x'  =  x,  x^  =  y,  and  x*  =  z.  Then,  for  instance,  Bj=  and  =  9B2/3x* 
=  dBy/dz.  We  further  employ  the  single  slash  notation  "  / "  to  denote  Euclidean  differentiation 
and  the  double  slash  "  // "  to  denote  non-Euclidean  differentiation,  where  the  curvature  of  the 
metric  space  must  be  considered.  For  the  moment,  we  will  be  concerned  only  with  Euclidean 
space,  for  which  our  example  reduces  to:  =  3By/3z  =  Bj^.  So,  in  general,  the  vector  and 

tensor  components  of  the  magnetic  field  are  denoted  as  B^  and  respectively,  where  the 
Greek  indices  |i  and  v  range  from  1  to  3.  As  a  further  notational  point,  we  emphasize  that  in  a 
Euclidean  or  Cartesian  coordinate  system,  there  is  no  distinction  between  contravariant  vector 
and  tensor  components  with  raised  indices  (e.g.,  )  and  covariant  vector  and  tensor 

components  with  lowered  indices  (e.g.,  )•  Consequently,  we  will  use  the  convention  that  all 

indices  are  to  remain  lowered  unless  they  are  being  summed  over.  Then,  in  accordance  with 
Einstein  summation  notation,  the  index  being  summed  over,  which  must  necessarily  appear 
twice  in  a  given  mathematical  expression,  will  have  one  index  raised  and  one  lowered.  Then,  as 
an  example: 

=  i  (17) 

This  convention  will  be  used  exclusively  unless  otherwise  explicitly  stated  to  the  contrary. 
Using  this  tensor  notation,  which  corresponds  to  that  of  Adler,  Bazin,  and  Schiffer  (1975), 


9 


eq.  (13)  becomes: 


—  —  C>/ji 

(18) 

while  eq.  (14)  becomes: 

< 

II 

CD 

< 

(19) 

Also,  recalling  that  the  point  of  observation  r  is  assumed  to  be  outside  of  the  magnetized  volume 

V  so  that  H(r)  =  B(r),  the  magnetostatic  eqs.  (4a)  and  (4b)  reduce  to  the  following: 

=0 

(20a) 

CD 

< 

1 

CD 

< 

II 

o 

(20b) 

Using  eq.  (19),  eqs.  (20a)  and  (20b)  may  also  be  written  as: 

=  0 

(21a) 

II 

(21b) 

Furthermore,  since  is  the  gradient  of  a  vector,  it  is  by  definition  a  tensor.  Equation  (21a) 

states  that  the  trace  (i.e.,  the  diagonal  element  sum)  of  the  tensor  is  zero,  while  eq.  (21b)  says 
that  this  tensor  is  also  symmetric.  These  two  symmetry  conditions,  which  may  be  viewed  as 
constraints  on  the  magnetic  field  gradient  tensor,  reduce  the  number  of  the  tensor's  independent 
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components  from  9  to  S.  Designers  of  magnetic  gradiometers  take  this  fact  into  account  to 
reduce  the  electronic  complexity  of  the  instrument.  However,  some  redundency  may  be  useful 
in  actual  field  situations. 

2.2  Related  Magnetic  Vectors,  Tensors,  and  Rotational  Invariants 

From  the  vector  and  the  gradient  tensor  ,  a  variety  of  potentially  useful  vectors, 

tensors,  and  scalar  invariants  can  be  constructed.  These  invariants  can  be  used  as  additional 
tools  to  help  identify  and  interpret  the  underlying  crustal  geology.  For  instance,  the  following 
vector  may  be  constructed; 


(22) 

The  vector  (units:  nTVcm)  may  be  viewed  as  the  projection  of  the  magnetic  gradient  tensor 
onto  the  magnetic  induction  vector.  From  it  may  be  constructed  its  characteristic 
rotational  invariant  (units:  nT^/cm): 

P  =  (23) 

Additionally,  limiting  the  discussion  only  to  tensors  of  rank  two,  it  is  possible  to  construct  the 
following  four  new  tensors: 

Piiv  3  P^pv  (24a) 

^iiv  3  p^  Bv  (24b) 
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Sx,M 


(24c) 


Bnv  ■  Bn  Bv  (24d) 

Similar  toisors  can  also  be  constructed  from  the  electric  and  gravitational  fields.  Three  of  these 
tensors,  (units:  nT^/cm^),  (units:  nTVcm^),  and  (units:  nT^)  are  symmetric,  while  the 
tensor  (units:  nTVcm)  has  both  symmetric  and  antisymmetric  parts.  By  taking  the  trace  of 

these  tensors,  several  additional  rotational  invariants  can  also  be  constructed.  Although  we  will 
not  dwell  upon  them,  we  point  them  out  so  that  these  invariants  may  be  exploited  for  their 
long*range  and  short-range  dependencies,  as  the  needs  of  a  particular  application  dictate.  We 
define  these  invariants  as  follows: 


p  = 

(25a) 

2  m 

(25b) 

SPm 

(25c) 

B  ■ 

(25d) 

The  last  invariant,  B,  is  just  the  total  magnetic  intensity,  while  SP  is  the  total  magnetic  gradient. 


From  these  invariants,  it  is  possible  to  generate  three  other  invariants  that  have  dimensions  of 
inverse  length  (i.e.,  units:  cm  '): 


Qi  ■  S^/B 


(26a) 


Q2  ■  p/B*  (26b) 

Qa  ■  (26c) 

It  so  happens  that  all  three  of  these  invariant  functions,  in  the  far  field,  vary  inversely  with 
distance  from  the  source  (i.e.,  as  f ').  This  can  be  seen  by  noting  that  in  the  far  field,  varies 
as  f  while  varies  as  f  ^ .  By  taking  ratios  of  these  invariants,  three  more  dimensionless 

invariants  can  be  generated  which  in  the  far  field  are  constant,  but  in  the  near  field  are  highly 
variable  and  may  be  quite  sensitive  to  subtle  variations  in  the  source  structure  or  composition. 
These  are; 


qj  «  Qi/Qj 

(27a) 

q2  *  Q2/Q3 

(27b) 

qs  ■  Q3/Q1 

(27c) 

Some  of  these  invariants  may  prove  to  be  more  useful  than  others  as  tools  for  geological 
interpretation.  One  that  does  seem  to  be  particularly  interesting  is  Q,,  which  is  simply  the  ratio 
of  the  total  magnetic  gradient  intensity  to  the  total  magnetic  field  intensity.  It  has  the  desirable 
feature  of  being  independent  of  the  orientation  of  the  magnetization  vector  of  the  source  body. 
Thus*  the  source  body  lateral  extention  is  more  sharply  defined  than  with  traditional  methods. 
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2  J  The  Uniformly  Magnetized  Rectangular  Prism 

The  static  magnetic  field  generated  by  an  arbitrarily  shaped  body  which  has  an  arbitrary 
distribution,  orientation,  and  intensity  of  magnetization  can  always  be  represented  as  having 
been  generated  by  a  collection  of  uniformly  magnetized  rectangular  prisms.  It  is,  of  course, 
possible  to  use  other  elementary  geometries  such  as  multifaceted  polyhedrons  as  the  basic 
building  blocks.  However,  for  local  field  analysis  we  see  no  special  advantage  in  doing  so. 
Instead,  we  prefer  to  maintain  the  unique  compatibility  that  rectangular  prisms  have  with  the 
Cartesian  coordinate  system. 

Now,  consider  a  single  rectangular  prism  with  sides  parallel  to  the  coordinate  axes,  the 
orientation  of  which  we  have  previously  defined.  The  center  of  this  prism  with  respect  to  the 
origin  of  the  coordinate  system  is  located  at  the  point  Tq  =  (x^,  yo,  z^).  The  prism  dimensions  are 
taken  to  be  X, ,  Xy ,  and  X^  along  their  respective  coordinate  axes.  When  this  prism  is  assigned  a 
uniform  magnetization  M,  having  constant  components  M,  =  M„  My  =  Mj ,  and  =  Mj  so 
that: 


M  =  M,l  +  Myj  +  Mzk  (28) 

then  it  is  possible,  as  is  shown  in  Appendix  A,  to  evaluate  eq.  (12)  for  the  scalar  potential,  0(r), 
in  closed  form  and  thereby  obtain  the  following  expression: 

0(r)=-Q^e“Mx  (30) 

where  the  elements  of  the  ft(r)  matrix,  which  are  explicitly  listed  in  Table  I,  are  purely 
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Table  1.  Elements  of  the  Matrix 


n„  =  .Kx-x').an-'[i2fgg^]  1.,,..,, 
ni2  =  -  (y  -  y')  ln[R  +  (z  -  z')]  1,/  .y-  ,*/ 

f2i3  =  -  (z  -  zO  ln[R  +  (y  -  yO]  L'.y'.z' 

^21  =  -  (X  -  XO  ln[R  +  (Z  -  Z')]  Ix'.y'.y' 

fl22  =  +  (y  - /)  L-.y,.,, 

fl23  =  -  (Z-z')ln[R  +  (x-x')]  L'.y'.x' 
flai  =  -  (X  -  x')  ln[R  +  (y  -  y')]  fx'.y'.z' 

^32  =  -  (y  -  yO  ln[R  +  (X  -  X')]  Ix'.y'.x' 

£2,3  =  +(2-z0tan-'[ii^^]  lx^y^x/ 
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geometric  in  character,  depending  only  on  the  position  r,  the  prism  dimensions,  and  the  pnsm 
location  .  Both  r  and  Tq  are  referenced  to  a  conveniently  chosen  origin. 

As  a  matter  of  notational  convenience,  we  have  defined  the  constant  vector  e: 

e  =  eii  +  £2]  +  e^k  (31) 

so  that  each  of  its  components  is  equal  to  one.  That  is; 

£i  =  £2  =  63  *  1  (32) 

The  virtue  of  defining  the  vector,  e,  in  this  manner  is  that  the  mathematical  forms  of  various 
geopotential  field  equations  are  considerably  simplified. 

In  Table  1,  the  notation  [  y.  is  used  to  indicate  that  the  primed  coordinates  are  to  be 
evaluated  at  the  boundaries  of  the  prism  as  indicated  in  Appendix  A.  Also  in  Table  l,the 
parameter  R  is  defined  to  be  the  Euclidean  distance  between  the  source  prism  and  the  point  of 
observation; 

R  =  y(x-xO^  +  (y-yO^  +  (z-z')^  (33) 

Although  it  is  possible  to  derive  an  expression  for  the  magnetic  induction  vector,  B(r),  by 


directly  taking  the  gradient  of  eq.  (30),  it  turns  out  to  be  simpler  to  combine  eqs.  (12)  and  (13) 
and  then  to  evaluate  the  resulting  integrals  as  shown  in  Appendix  B.  The  vector  components  of 
the  magnetic  induction  vector,  thus  obtained,  are  expressed  quite  simply  in  tensor  notation  as; 


Bh  =  A/ Ml 


(34) 


whore  the  elements  of  the  3  x  3  matrix  A(r),  which  are  listed  in  Table  2,  are,  like  those  of  the 
matrix  Q  to  which  it  is  related,  purely  geometrical  in  character. 

Figures  (la),  (lb),  and  (Ic)  illustrate  the  three  components  (x,  y,  and  z  respectively)  of  the 
magnetic  field  vector  B  generated  according  to  eq.  (34)  by  a  single  prism  with  unit 
magnetization  M  oriented  vertically  downward  in  the  positive  Z  direction.  The  computations 
were  performed  on  a  surface  situated  a  unit  distance  above  the  prism.  Figure  (Id)  illustrates  the 
magnetic  total  intensity  for  the  same  situation.  The  scale  in  these  figures  is  arbitrary.  Therefore, 
these  figures  could  represent  the  magnetic  field  generated  by  a  rock  sample  a  few  centimeters  in 
diameter  in  the  laboratory,  or  they  could  i«present  the  magnetic  field  generated  by  a  large 
volume  of  ocean  crust  several  kilometers  in  diameter.  Figures  (2a),  (2b),  and  (2c)  are  vector 
compments  of  the  magnetic  induction  B,  (North),  By  (East),  and  B,  (Vertically  Down) 
req)ectively  for  a  multiple  prism  situation  simulating  a  region  encompassing  a  magnetic  reversal 
in  the  oceanic  crust.  The  magnetic  reversal  axis  straddles  the  East-West  coordinate  axis.  North 
of  this  axis,  the  magnetization  is  vertically  upward,  while  South  of  this  axis,  the  magnetization  is 
vertically  downward.  The  magnetizations  are  of  unit  magnitude  in  both  directions.  The  prism 
dimoisions  are  arbitrarily  chosen  to  be  of  unit  length.  Figure  (2d)  is  the  total  magnetic  intensity 
for  this  situation.  The  observation  plane  is  a  unit  distance  above  the  line  of  prisms.  The 
magnetic  field  contribution  of  each  prism  is  added  vectorially  at  each  grid  site  on  the 
observational  plane. 
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Figure  la.  North-X  Component  of  a  Vertically  Magnetized  Prism, 
(units:  nT,  arbitrary  scale) 


Figure  lb.  F-^st-Y  Component  of  a  Vertically  Magnetized  Prism, 
(units:  nT,  arbitrary  scale) 
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Figure  Ic.  Vertical-Z  Component  of  a  Vertically  Magnetized  Prism, 
(units:  nT,  arbitrary  scale) 
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Figure  Id.  Total  Intensity  of  a  Vertically  Magnetized  Prism, 
(units:  nT,  arbitrary  scale) 


Figure  2a.  North-X  Component  in  a  Region  of  Vertical  Magnetization  Reversal, 
(units:  nT,  arbitrary  scale) 
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Figure  2b.  East-Y  Component  in  a  Region  of  Vertical  Magnetization  Reversal, 
(units:  nT,  arbitrary  scale) 


Figure  2c.  Vertical-Z  Component  in  a  Region  of  Vertical  Magnetization  Reversal, 
(units:  nT,  arbitrary  scale) 


»  .• 

I  !  • 
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Figure  2d.  Total  Intensity  in  a  Region  of  Vertical  Magnetization  Reversal, 
(units:  nT,  arbitrary  scale) 


It  is  perhaps  helpful  to  view  the  12  matrix  as  a  transformation  that  maps  the  magnetization 
vector  M  into  the  magnetic  potential  <I),  while  the  matrix  A  acts  as  a  transformation  that  maps 
the  magnetization  vector  into  the  magnetic  induction  B.  The  relationships  among  components 
of  these  two  matrices,  which  are  in  fact  tensors  of  rank  two,  are  found  by  combining  eqs.  (18), 
(29),  and  (34),  which  yields: 

A(iv  —  linff/vE*  (35) 

The  magnetic  gradient  due  to  the  rectangular  prism  is  obtained  by  combining  eqs.  (19)  and 
(34),  which  yields: 

=  A/^Mx  (36) 

The  gradient  of  A  appearing  in  eq.  (36)  forms  a  tensor  of  rank  3.  It  maps  the  magnetization  of 
the  rectangular  prism  into  the  magnetic  gradient  tensor.  The  components  of  this  third-rank 
tensor  are  listed  in  Table  3.  These  components  are  obtained  by  taking  the  indicated  derivatives 
of  the  elements  of  the  A  matrix  and  simplifying  the  resulting  algebraic  forms. 

The  A  matrix  is  symmetric,  and  if  the  point  of  observation  r  extends  beyond  the  boundary  of 
the  volume  V  (denoted  as  3V  ),  then  it  is  also  traceless.  If  r  is  within  this  boimdary,  then  the 
trace  of  A  is  a  nonzero  constant.  These  results  may  be  summarized  in  the  following  way: 


A|iv  =  A 


VM 


(37a) 


I  -An  r  <  3V' 
I  0  r  >  av' 


(37b) 
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Wc  are  primarily  concerned  with  the  situation  for  which  the  trace  is  zero.  Together,  these  last 
two  equations  imply  that  A  has  just  5  independent  components.  Using  eqs.  (37a)  and  (37b),  it  is 
possible  to  deduce,  in  conjunction  with  Table  3,  the  following  symmetry  relations  among  the 
gradient  components  of  A: 


A^ylX  —  Avll/X 

(38a) 

A)jv/x  —  Axii/v  —  Avx/n 

(38b) 

A'.a  -  0 

(38c) 

These  symmetry  relations  substantially  reduce  the  number  of  independent  components  of  the 
gradient  of  the  A  tensor  from  a  possible  27  to  just  7. 

2.4  Geomagnetic  Inversion 

Still  considering  a  single  prism,  it  should  be  clear  that  at  any  observation  point  r,  eq.  (34)  can 
be  formally  inverted  to  give: 

Mx  =  Ax'"  (39) 

where  A  is  the  inverse  of  the  3x3  matrix  A  and  has  elements  that  satisfy  the  following  relation: 
a/Aiv  =  5,v  (40) 

where  is  the  Kronecker  delta  function.  Subsequent  insertion  of  eq.  (39)  into  eq.  (36)  yields  a 
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relationship  between  the  magnetic  field  and  its  gradient  that  is  independent  of  the  magnetization: 


(41) 

This  is  the  main  result.  Essentially,  a  decoupling  of  the  geophysical  and  geometric 
parameters  has  taken  place.  Since  the  magnetic  field  and  its  gradient  are  both  measurable 
quantities,  eq.  (41)  is  in  essence  a  nonlinear  system  of  five  independent  equations  for  the  six 
unknown  geometric  parameters  of  the  magnetized  prism,  namely:  and  X, ,  ,  and  Since 

there  are  more  unknowns  than  equations,  this  system  of  equations  is  underdetermined.  For  this 
reason  and  also  because  magnetic  field  measurements  contain  errors  and  extraneous 
environmental  noise,  it  is  necessary  to  make  some  additional  assumptions  regarding  the  size  or 
location  of  the  prism  in  order  to  reduce  the  number  of  unknown  geometric  parameters. 
Alternatively,  it  is  possible  to  make  many  field  measurements  at  a  variety  of  locations,  in 
which  case,  this  system  of  equations  will  be  overdetermined.  In  either  case,  or  when  the  two 
approaches  are  combined,  the  problem  is  a  nonlinear  one  which  can  be  solved  via  stochastic 
inversion  techniques.  Having  thus  determined  the  prism's  geometric  parameters,  the  tensor  A(r) 
is  considered  to  be  known  and  is  unique  in  the  least-squares  sense  as  being  generated  by  the 
"most  likely"  or  "optimal"  set  of  prism  parameters.  Knowledge  of  this  tensor  then  permits  eq. 
(39)  to  be  solved  via  a  straightforward  linear  inversion  to  obtain  the  "most  likely"  or  "optimal" 
magnetization  (the  geophysical  parameter)  of  the  prism,  which  is  referred  to  as  the  "equivalent" 
magnetization  of  the  prism.  This  is  the  magnetization  one  expects  to  obtain  if  the  actual 
magnetization  were  truly  smeared  out  uniformly  over  the  entire  prism,  and  the  "most  likely" 
geometric  parameter  solution  is  also  the  environmentally  "true"  solution. 
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2.5  PractioU  Applications  to  Large  Area  Aeromagnetic  Surveys 

To  determine  the  crustal  magnetization  and  depth  to  the  magnetic  basement  (i.e.,  to  the  basalt 
layer),  it  is  not  actually  necessary  to  measure  the  gradient  Held  in  order  to  apply  the  inversion 
technique  described  above.  If  a  2-dimensional  vector  aeromagnetic  survey  is  performed  over  a 
given  area,  it  is  possible  to  use  the  vector  magnetic  data  to  compute  a  rectangular  harmonic 
potoitial  function  from  which  the  gradient  data  may  be  computed.  The  procedures  used  in  the 
July  1981  Juan  de  Fuca  Project  MAGNET  vector  aeromagnetic  survey  performed  by  the  Naval 
Oceanognq)hic  Office  serve  as  a  practical  example  of  the  kind  of  data  manipulation  that  is 
necessary  to  use  this  inverse  modeling  technique  successfully  (Quinn  and  Shiel  [1993]). 

The  Juan  de  Fuca  survey  area,  which  covered  the  geographical  region  from  AT  N  latitude  to 
51^  N  latitude  and  from  124**  W  longitude  to  130‘’  W  longitude,  was  densely  surveyed  at  a  SOO-ft 
altitude  in  the  East- West  direction  with  approximately  3  nautical  miles  between  these  survey 
tracks.  A  few  North-South  tracks  were  also  flown  for  control.  The  effective  along-track  data 
sanq)le  rate  was  0.5  Hertz,  while  the  speed  of  Project  MAGNETs  RP-3D  Orion  aircraft  was 
^proximately  240  nautical  miles  per  hour.  Magnetograms  from  a  nearby  geomagnetic 
observatory  at  Victoria,  British  Columbia,  were  used  to  monitor  and  remove  temporal  magnetic 
variations.  A  portable  Vector  Magnetic  Ground  Station  (VMGS)  was  also  established  at 
Me  Chord  Air  Force  Base  near  Tacoma,  Washington,  for  the  same  purpose.  This  data  set  was 
edited  for  spurious  elev'tronic  noise  spikes,  and  the  Main  magnetic  field  was  removed  using  the 
1980  Epoch  Goddard  Space  Flight  Center  (GSFC  12/83)  spherical  harmonic  model  (Langel  and 
Estes  [1985])  up  to  harmonic  degree  12.  The  residual  X-,  Y-,  and  Z-components  of  the  entire 
data  set  were  each  uniformly  gridded  and  interpolated  as  necessary  into  lat  x  Ion  cells  of  1.5 
areminutes  x  1.5  areminutes  (i.e.,  roughly  half  the  East- West  track  line  separation)  such  that 
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each  magnetic  component  consisted  of  2S6  x  2S6  grid  points.  The  gridding  process  has  the 
advantage  of  smoothing  out  instrument  noise,  navigational  errors,  and  other  environmental  noise 
that  cannot  otherwise  be  accounted  for  and  thus  eliminated.  The  coordinates  of  each  grid  point 
were  then  transformed  from  geodetic  to  rectangular  coordinates,  using  the  spherical  coordinate 
system  as  an  intermediate  step.  At  the  same  time,  the  vector  components  at  each  grid  point  were 
rotated  from  geodetic  to  rectangular  coordinates,  also  using  the  spherical  coordinate  system  as  an 
intermediate  stage.  Biases  and  linear  trends  in  the  North-South  and  East-West  directions  were 
then  determined  and  removed  from  the  residual  grids  of  their  corresponding  components.  These 
biases  and  trends  were  also  saved  as  the  "regional"  models  of  the  survey  area.  Subsequently,  a 
2-dimensional  Fast  Fourier  Transform  (2-D  FFT)  was  performed  on  only  the  Z-component  grid. 
This  procedure  results  in  a  set  of  complex  Fourier  coefficients  in  the  waveniimber  domain  which 
are  algebraically  related  in  a  simple  way  to  the  real  coefficients  of  a  rectangular-harmonic 
magnetic  potential-function  model  of  the  surveyed  area.  This  model  is  composed  of  2S6  x  256  = 
65,536  coefficients.  Having  determined  the  rectangular  harmonic  model  in  this  way,  it  is  a 
straightforward  process  to  recompute,  using  the  2-D  FFl',  not  only  the  vector  magnetic 
components  but  the  magnetic  gradient  components  as  well.  The  result  is  a  set  of  3  vector 
component  grids  and  5  independent  gradient  component  grids,  each  with  256  x  256  grid  points. 
Since  all  of  these  magnetic  component  grids  were  derived  from  the  original  Z-component  grid, 
they  are  necessarily  consistent  with  each  other.  The  original  vector  component  grids  derived 
directly  from  the  survey  data  would  not  have  been  as  consistent  because  each  vector  component 
obsoration  contains  unique  instrumental  errors  and  environmental  noise  which  are  then  passed 
on  to  the  grids  generated  from  those  measurements.  Therefore,  in  the  inversion  process  that 
follows,  we  use  the  X-  ,  Y-,  and  Z-component  grids  generated  from  the  model,  which  in  turn  is 
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derived  from  just  the  observed  Z-component  gridded  data,  rather  than  using  the  X-  and 
Y-component  grids  generated  directly  from  the  observed  data.  The  computed  and  observed 
Z-component  grids  are  identical  since  this  component  is  used  to  generate  the  model. 

At  this  stage  each  grid  point  has  associated  with  it  the  "observed"  magnetic  vector 
components  and  the  "observed"  magnetic  gradient  components  .  Here,  the  word 

"observed"  is  used  loosely  since  each  grid  value  of  a  magnetic  vector  or  gradient  component  is 
based  indirectly  on  measured  data  through  the  rectangular  harmonic  model,  which  in  turn  is 
derived  from  observed  data.  Magnetic  gradients,  generated  from  eq.  (41)  using  "only"  the 

"observed"  magnetic  "vector"  components  will  be  referred  to  as  the  "computed"  magnetic 
gradients.  Thus,  for  each  magnetic  gradient  component  at  each  grid  point  there  exists  an 
"observed"  and  a  "computed"  value.  It  is  therefore  possible  to  solve  for  the  geometric 
parameters  for  each  grid  point  by  setting  up  a  function  which  is  the  sum  of  the  square  of  the 
differences  between  the  "computed"  and  "observed"  values  of  the  five  independent  gradient 
components  for  each  grid  point.  In  the  Juan  de  Fuca  case,  we  first  made  a  few  assumptions 
based  on  knowledge  of  magnetic  sources  in  the  oceanic  crust  to  constrain  the  inversion  problem 
in  order  to  eliminate  some  variables.  Since  it  is  known  that  the  oceanic  crust  typically  has  two 
magnetic  layers,  one  about  .5  km  thick  with  reasonably  high  magnetization  and  another  layer 
about  1  km  thick  directly  below  the  first  layer,  but  considerably  less  magnetic,  we  fixed  the 
thickness  of  a  single  prism  to  be  1 .5  km  thick.  The  lateral  prism  dimensions  were  also  fixed  to 
be  5  km  on  a  side  (about  the  same  length  as  the  estimated  depth  to  source).  The  prism  was 
centered  below  the  particular  grid  point  of  interest  so  that  there  was  a  minimum  of  2.5  km  from 
the  grid  point  to  the  edge  of  the  prism.  From  a  geomorphological  point  of  view,  dramatic 
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variations  in  depth  or  magnetization  should  not  dominate  for  such  small  lateral  distances. 
Consequently,  the  influence  of  magnetized  material  much  beyond  a  prism  edge,  1.5  km  from 
the  observation  point,  should  be  comparatively  small  and  can  therefore  be  neglected. 
Consequently,  the  only  parameter  yet  to  be  determined  is  the  magnetic  source  depth  z^  of  the 
prism's  center.  It  can  be  evaluated  by  minimizing  the  function: 

3  3 

~  ~  ^|lv  (oJw)  (^2) 

|l>l  v«l 

This  function  was  minimized  by  varying  the  depth  parameter  at  .1-km  increments  from  the 
known  bathymetry  base  on  NAVOCEANO's  5-Minute  Digital  Bathymetric  Data  Base  (DBDB5) 
downward,  well  beyond  the  poorly  estimated  Curie  depths,  computing  the  x^  function  at  each 
depth  and  noting  where  the  minimum  occurred.  Having  obtained  the  magnetic  source  depth,  , 
in  this  marmer,  eq.  (39)  was  used  to  compute  the  equivalent  magnetization  of  the  prism.  This 
procedure  was  performed  independently  for  each  individual  grid  point  until  the  entire  survey 
area  was  covered,  yielding  a  256  x  256  grid  for  magnetic  source  depth  and  the  same  size  grid  for 
each  of  the  three  components  of  the  magnetization  vector.  The  results  from  one  grid  point  to  the 
next  were  smoothly  varying  (i.e.,  not  noisy)  over  the  entire  area.  This  is  attributed  to  smoothing 
that  takes  place  during  the  gridding  process  and  also  to  the  additional  use  of  self-consistent 
magnetic  parameters,  all  derived  from  the  same  rectangular  harmonic  potential.  Since  the 
inverted  data  are  generated  on  the  same  grid  as  the  input  data,  the  results  are  easily  profiled  as 
functions  of  either  latitude  or  longitude. 
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Because  the  dimensions  of  the  prism  were  somewhat  arbitrarily  chosen,  the  resulting  depths, 
though  reasonable,  were  not  taken  as  absolute.  This  is  not  a  limitation  of  the  method,  since  we 
could  have  optimized  the  prism  dimensions  along  with  the  depth.  Prism  dimension  constraints 
were  applied  to  the  Juan  de  Fuca  survey  due  to  the  large  size  of  the  survey  area  and  the  resulting 
computational  burden  it  entailed  on  an  older  mainframe  computer.  Also,  some  additional 
non-uniqueness  concerns  arise  when  all  geometric  parameters  are  allowed  to  float  freely. 

Finally,  improved  depth-to-source  determinations  are  possible  through  more  detailed 
modeling  techniques  using  multiple  prisms  stacked  vertically  below  a  grid  point  and  which  are 
assumed  to  have  varying  magnetization,  by  iterating  (i.e.,  by  computing  the  residual  magnetic 
field  and  gradient  field  components  after  the  characteristics  of  one  prism  have  been  determined 
and  solving  for  the  depth  and  magnetization  of  a  second,  less-thick  prism  constrained  to  be  at 
shallower  depths  below  the  same  grid  point)  and  by  making  other  fairly  simple  refinements. 

2.6  Inverse  Magnetic  Modeling  with  Multiple  Prisms 

Although  reasonable  results  can  be  obtained,  even  for  large  ocean  areas  of  complicated 
geomorphology  such  as  the  Juan  de  Fuca  region,  improved  results  ought  to  be  expected  through 
more  detailed  modeling  with  multiple  prisms.  In  this  case  we  consider  a  collection  of 
magnetized  prisms  foliated  within  part,  or  all,  of  the  oceanic  crust  beneath  a  surveyed  area. 
Depending  on  the  situation,  one  may  fix  the  number  of  prisms  a  priori  or  view  the  number  of 
prisms,  N,  as  another  parameter  to  be  determined  as  part  of  the  least-square  or  stochastic 
minimization  procedure.  One  can  further  allow  the  prisms  to  vary  in  size  and  to  overlap,  in 
which  case  the  resulting  magnetizations  in  the  overlapping  regions  would  simply  add  vectorially. 
Alternatively,  one  could  constrain  the  prisms  not  to  overlap  and  perhaps  even  fill  the  entire 
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volume  of  the  oceanic  crust  below  the  siurveyed  area  with  adjacent  prisms,  each  of  the  same 
fbced  dimensions.  Regardless  of  the  choices  made,  the  multiprism  inversion  technique  is 
basically  the  same,  although  computationally  more  demanding,  than  that  for  a  single  prism. 

The  generalization  to  the  multiprism  case  is  straightforward.  It  is  now  assumed  that  there  is  a 
total  of  K  observations  of  the  magnetic  field  and  gradient  field  components  which  are  presumed 
to  be  the  composite  of  all  fields  generated  by  a  total  of  N  prisms.  Then,  at  the  k'th  observation 
point  ,  eqs.  (34)  and  (36)  generalize  to  the  following: 


«Hvk=  SfV/vl  Mx»  k=I.2....K  (44) 

a. I  ^  ''•‘O 

where  the  notation  has  been  chosen  to  convey  the  following  meanings:  *  B^(r J, 

■  ^^yC^k),  (\^)to*  \ Vk.  f'n).  )kn  *  ly  (^.O-  Additionally,  is  the 

X-component  of  magnetization  of  the  n'th  prism,  which  is  centered  about  the  point  r'o„. 
Equations  (43)  and  (44)  can  be  written  more  succinctly  in  matrix  form  as  follows: 

J»R=T/inx  (45) 

=  'P/vilIx  (46) 

where,  for  a  total  of  K  observations  and  N  prisms,  the  matrices  in  the  above  equations  take  the 
following  forms: 
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niji 


'  Mi, 
Mu 
Mu 


(51) 


I,  MiN  ) 


A  formal  inversion  of  eq.  (45)  then  gives  the  magnetization  matrix  iJli  in  terms  of  the  matrix 
containing  the  magnetic  field  vector  component  observations: 

i«x  =  nivr''‘‘)B„  (52) 

where  y'"'*  is  the  transpose  of  T''**  with  respect  to  the  Latin  indices  n  and  k,  and  where  the  matrix 
Hxv  is  the  inverse  of  the  N  x  N  matrix  rivXt  which  is  defined  as  follows: 

Ovx  =  fv**  (53) 

Inserting  eq.  (52)  into  eq.  (46)  accomplishes  the  desired  separation  of  the  geophysical  and 
geometric  parameters  and  yields: 


=  'P/vnx„f“‘*i5p  (54) 

The  right  side  of  this  equation,  by  virtue  of  the  three  Einstein  summations  involving  the  indices 
a,  P,  and  K  contains  27  terms,  each  term  being  the  product  of  4  matrices.  Equation  (54)  is 
therefore  a  rather  large  system  of  nonlinear  equations  involving  only  the  geometric  parameters 
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associated  with  the  N  prisms  and  must  be  solved  iteratively  using  stochastic  or  generalized 
inversion  techitiques.  Thereby,  a  unique  solution  in  the  "optimized"  least-squares  sense  can  be 
obtained.  The  number  of  prisms,  N,  can  also  be  included  as  one  of  the  unknown  parameters  or 
die  number  of  prisms  may  be  set  a  priori.  The  solution  for  the  geometric  parameters,  thus 
determined  from  eq.  (S4),  then  provides,  in  conjunction  with  eq.  (52),  a  straightforward  linear 
inversion  problem  for  the  determination  of  the  vector  components  of  the  magnetization  in  each 
of  the  N  prisms. 

2.7  The  Computational  Algorithm  MAGREP 

Subroutine  MAGREP  (MAGnetic  field  due  to  a  REctangular  Prism)  computes  the  magnetic 
scalar  potential  d>(r)  from  eq.  (30),  the  magnetic  induction  vector  B(r)  from  eq.  (34),  and  the 
nine  elements  of  the  magnetic  gradient  tensor  ^r)  from  eq.  (36)  given  the  magnetization, 

position,  and  dimensions  of  prism.  It  is  assumed  that  the  user  of  this  routine  has  written  a  driver 
program  for  the  MAGREP  subroutine  that  defines  the  prism  parameters  and  passes  them  to 
MAGREP  through  a  common  block  called  /MAGBLK/.  It  is  further  assumed  that  the  position 
of  the  prism  is  referenced  to  an  origin  located  at  the  lower  left-hand  comer  of  a  user-defined 
surface.  The  point  of  observation  is  also  referenced  to  this  origin.  The  rectangular  components 
of  the  observation  point  are  passed  to  MAGREP  through  the  CALL  statement,  which  also 
returns  the  computed  magnetic  field  values.  The  field  observed  at  a  single  observation  point 
genmted  by  multiple  prisms  is  the  sum  of  outputs  of  MAGREP  for  each  prism.  The  prism  can 
be  rotated  with  respect  to  the  rectangular  coordinate  axes  associated  with  the  coordinate  origin 
through  a  set  of  three  Euler  angles  that  are  also  passed  to  MAGREP  from  the  driving  program 
through  the  common  block.  When  the  Euler  angles  are  zero,  the  prism  is  unrotated,  and  its  sides 


are  parallel  to  the  coordinate  axes.  Documentation  with  a  detailed  description  of  these  and  other 
parameters  is  contained  in  the  FORTRAN  code  of  the  subroutine  which  is  listed  in  Appendix  E. 
This  subroutine  is  well  suited  for  simulating  magnetic  fields  generated  by  a  variety  of 
geophysical  entities  such  as  seamounts,  as  well  as  those  of  man-made  origin,  such  as  ships, 
planes,  trains,  etc. 

3.  THE  GRAVITY  FIELD 
3.1  Theoretical  Background 

Modeling  a  gravitational  field  that  is  generated  by  some  crustal  feature  or  collection  of 
features  begins  with  a  pair  of  equations  that  may  be  considered  as  the  gravitational  analogue  of 
Maxwell's  magnetostatic  equations.  The  static  gravitational  field  equations  are: 


V*g  s=  -4xGp 

(55a) 

Vxg  *  0 

(55b) 

vdioe  g  is  the  gravitational  acceleration  vector  (units:  gal  =  cm/sec^ ),  G  is  the  universal 
gravitational  constant  (G  =  6.6720  x  10  *  cmVgm-sec^),  and  p  is  the  density  contrast  of  the 
Earth's  crust  (units:  gm/cm* ). 

Equation  (SSb)  indicates  that  the  gravitational  acceleration  vector  may  be  considered  as  the 
gradient  of  some  scalar  potential  U(r}  (units:  cmVsec^ ),  the  units  of  which  correspond  to  those 
of  energy  per  unit  mass.  Consequently: 

g=-VU  (56) 
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This  equation,  when  combined  with  eq.  (SSa),  yields  the  familiar  Poisson  equation  for  the 
gravitational  field: 

V*U  «  4xGp  (57) 

>i^di  has  the  well-known  solution; 

tSt*'’'"  (5«) 

Equation  (58)  is  valid  for  any  observation  point  r  regardless  of  whether  or  not  r  is  inside  or 
outside  the  gravity  field  source  volume  V.  In  the  event  that  the  observation  point  r  is 
outside  of  this  volume,  then  eq.  (58)  must  also  satisfy  Laplace's  equation; 

V^U  =  0  (59) 

The  solutions  to  this  equation  are  also  well-known,  depending  on  the  coordinate  geometry  as 
rectangular  harmonic  functions,  spherical  harmonic  functions,  etc. 

The  gradient  of  the  gravity  acceleration  vector,  g,  yields  the  gravity  gradient  tensor 
(units:  sec'^): 

»(r)  =  Vg(r)  (60) 

Both  g(r)  and  ^r)  are  measurable  quantities.  The  state  of  the  art  in  gravity  gradiometry  is 
discussed  by  Jordan  (1978,  1985)  and  Wells  (1983). 
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As  in  the  magnetic  case,  the  gravity  acceleration  vector  and  the  gravity  gradient  tensor  may 
be  decomposed  into  their  respective  Cartesian  components  such  that: 

g(r)  =  g»(r)i  +  gy(r)i  +  gz(r)k  (61) 

and 


»(r)  = 


f  ill 

agi 

dx 

dx 

dx 

9g« 

3gy 

dgz 

9y 

ay 

3gx 

agy 

dgz 

^  9z 

dz 

dz 

(62) 


In  tensor  notation,  the  components  of  the  gravity  acceleration  vector  are  denoted  as  g^,,  while 
the  components  of  the  gravity  gradient  tensor  are  denoted  as  ,  where,  as  in  the 
magnetic  case,  the  Greek  subscripts  p,  v,  X,  etc.,  range  from  1  to  3  so  that  gi  =  g, ,  82  “  gy . 
gj  =  g,  ,  and  ^f2  =  3g*/9y,  etc.  Using  this  notation,  eqs.  (55a)  and  (55b)  are  cast  into  the 
following  forms: 

g**/^  =  -4nGp  (63a) 

git/v  -  gv/|»  =  0  (63b) 

Also,  eq.  (56)  becomes: 

gii  =  -U/n  (64) 
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eq.  (S7)  takes  the  form: 


=  4rGp 


(65) 


Finally,  eq.  (60)  in  tensor  notation  becomes: 

=  Sit/v  (66) 

Now,  by  combining  eq.  (66)  with  eqs.  (63a)  and  (63b)  it  is  found  that: 

=  -4nGp  (67a) 

=  «v,i  (67b) 

That  is,  the  gravity  gradient  tensor  is  symmetric,  and  its  trace  is  proportional  to  the  density  p  at 

the  point  of  observation  r.  In  particular,  when  r  is  outside  the  volume  V*,  the  trace  is  zero,  so 
that: 

=  0  (68) 

Finally,  we  know  that is  a  tensor  because  it  is  the  gradient  of  the  gravity  acceleration  vector, 

and  by  definition,  the  gradient  of  any  vector  or  tensor  is  also  a  tensor. 
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3.2  Rdttcd  Gravity  Vectors,  Tensors,  and  Rotational  Invariants 


As  in  the  magnetic  case,  the  components  of  the  gravitational  acceleration  vector  and  the 
components  of  the  gravity  gradient  tensor  may  be  combined  in  various  ways  to  form  a 

variety  of  other  vectors,  tensors,  and  rotational  invariants,  in  direct  analogy  with  the  magnetic 
field.  These  quantities  are  additional  tools  that  may  be  used  to  examine  the  field  observations  in 
ways  that  hitherto  have  been  overlooked  or  which  may  become  more  useful  as  sophisticated 
gradiometric  instruments  become  more  widely  available.  They  can  yield  supplementary 
information  concerning  the  underlying  structure  and  density  contrast  of  a  surveyed  area  that  may 
be  used  to  constrain  the  density  morphology  model  of  that  area.  From  the  vector  g^,  and  the 
tensor  we  can,  for  instance,  construct  the  vector  \  and  its  corresponding  scalar  invariant  y 
(units:  cm/sec^ ): 


TTm  ■ 

(69a) 

Y-  77^ 

(69b) 

Additionally,  the  following  four  tensors  of  rank  two  may  also  be  constructed: 


Yitv  ■  YhYv 

(70a) 

■  Tr  gv 

(70b) 

^RV  * 

(70c) 

gRv  ■  gRgv 

(70d) 
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Three  of  these  tensors,  (units:  cm^/sec*  ).  (units:  sec  “*  ),  and  (units:  cmVsec^  )  are 
symmetric,  while  the  tensor  (units:  cm^/sec* )  has  both  symmetric  and  antisymmetric  parts. 
Rotational  invariants  constructed  from  the  traces  of  these  tensors  are  defined  as  follows: 


7" 

(71a) 

^  m 

(71b) 

9^  *  ii 

(71c) 

g 

(71d) 

The  last  invariant,  g  (units:  cm/sec^ ),  is  simply  the  total  gravitational  acceleration  (i.e.,  the 
magnitude  of  the  gravitational  acceleration  vector),  while  ^  (units:  sec'^  )  is  the  total 
gravitational  gradient.  Equations  (69b)  and  (71a)  are  different  expressions  for  the  same 
parameter,  y.  The  invariant  ^  has  the  same  units  as  the  tensor  itself 

In  the  far  field,  the  gravitational  acceleration,  g,  will  vary  as  l/R^  where  R  is  the  distance 
from  the  source  of  the  gravitational  field,  while  the  gravitational  gradient  ^  will  vary  as  1/R* . 

Consequently,  in  the  far  field,  y  will  vary  as  1/R* ,  and  ^  will  vary  as  1/R^  .  Furthermore, 

several  combinations  of  these  invariant  parameters  form  other  invariant  parameters  that  have  a 
1/R  dependence  in  the  far  field  and  also  have  units  of  inverse  length  (i.e.,  units  of  cm  ' ).  These 
are: 
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J 


Pi  5  ^/g 


(72a) 


?2  3  ^/g^ 

(72b) 

III 

0U 

(72c) 

Ratios  of  these  parameters  in  turn  will  yield  dimensionless  parameters  that  are  constants  in  the 

far  field,  but  which  may  be  quite  sensiti'^e  to  subtle  variations  in  the  source  geometry  or 

composition  in  the  near  field.  These  are: 

Pi  =  P1/P2 

(73a) 

P2  =  P2/P3 

(73b) 

P3  =  P3/P1 

(73c) 

These  parameters  simply  provide  a  different  way  of  looking  at  the  gravitational  field  and  ought 
to  be  viewed  as  additional  tools  in  the  repertoire  of  modem  geopotential  analysis,  which  may, 
under  the  appropriate  circumstances,  be  exploited.  Their  value  stems  from  the  fact  that  their 
dependence  on  the  radial  distance  from  the  source  is  quite  different  from  that  of  the  gravity  field 
or  its  corresponding  gradient  tensor  field. 

3.3  The  Uniformly  Dense  Rectangular  Prism 

The  last  section  reviewed  the  essence  of  Newtonian  gravitational  field  theory  without 
specifying  the  nature  of  the  gravitating  body  or  bodies.  Now,  as  in  the  magnetic  case,  we 


55 


rq»tsent  sources  of  gravitational  fields  as  being  composed  of  a  collection  of  rectangular  prisms, 
each  of  which  is  assigned  a  uniform  density,  which  may  be  different  from  one  prism  to  another. 

Considering  just  one  such  prism  oriented  with  its  sides  parallel  to  the  coordinate  axes,  we 
again  characterize  the  prism  as  having  its  center  at  the  point  =  (Xo ,  Vo .  Zq)  '’'^^h  respect  to 
some  coordinate  origin  and  as  having  dimensions  ,  Xy  ,  and  \  along  their  respective 
coordinate  axes.  Then,  given  that  the  density,  p,  of  the  prism  is  constant,  eq.  (S8)  can  be 
evaluated  for  the  gravity  potential  U(r),  as  is  shown  in  Appendix  C,  yielding  the  result; 


U(r)  =  -e»roxD>-  (74) 

where  the  elements  of  the  F  matrix  are  listed  in  Table  4,  where  represents  the  components 
of  the  constant  vector  e,  which  was  previously  defined  in  eqs.  (31)  and  (32),  and  where 
represents  the  components  of  the  density  vector  D  (units;  sec'^  ),  which  is  defined  by  analogy  to 
the  magnetization  vector  M  so  that; 

D  s  Gpc  (75) 

Then,  in  terms  of  its  components  ,  we  have; 

Djt  =  Gpex  (76) 

This  equation  simply  means  that  D,  =  =  D,  =  Gp  .  Introduction  of  the  density  vector  allows 

us  to  put  eq.  (74)  and  related  equations  into  a  convenient  and  compact  form  that  is  comparable 
to  eq.  (30)  for  the  magnetic  case.  The  F  matrix  (units;  cm^ )  is  symmetric,  as  can  be  directly 
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Table  4.  Elements  of  the  F  Matrix 


r„  =  +(x-xVtan-[^:g:^] 

r,2  =  -  (X  -  x')(y  -  y')  ln[R  +  (z  -  z')]  I »/ 

Fis  =  -  (x-xO(z-z')in[R  +  (y-/)]  L'.y'.*/ 
Fai  =  -  (y-yO(x-xOln[R  +  (z-z')]  L/.y/.x' 

r23  =  -  (y  -  yO{z  -  z')  ln[R  +  (x  -  x')]  I  K>,y',z' 
Fsi  =  -  (z  -  z')(x  -  x')  ln[R  +  (y  -  yO]  Ix'.y'.i' 
r32  =  -  (z  -  zO(y  -  y')  ln[R  +  (x  -  x')]  L'.y'.z' 
r33  =  +(z-z')^tan-'[ii^^]  L,.y..,, 
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verified  by  an  examination  of  Table  4.  Therefore: 

r^v  =  Fvii  (77) 

However,  it  has  a  nonzero  trace  regardless  of  whether  or  not  the  point  of  observation  is  inside  or 
outside  the  gravitational  source  volume  V  . 

Taking  the  gradient  of  eq.  (70)  will  yield  the  gravity  acceleration  vector  g  due  to  a  prism  of 
uniform  density.  It  turns  out  to  be  somewhat  easier  to  obtain  the  mathematical  expression  of  the 
acceleration  vector  by  taking  the  gradient  of  eq.  (S8)  and  subsequently  evaluating  the  resulting 
integrals,  as  is  shown  in  Appendix  D.  These  integrals  are  the  same  as  those  previously 
encountered  for  the  magnetic  potential.  Performing  the  necessary  computations  yields  the 
following  result  for  the  gravity  acceleration  vector  components: 

(78) 

where  the  elements  of  the  Q  matrix  (units:  cm)  have  been  previously  encountered  and  are  listed 
in  Table  1 .  This  matrix  has  no  special  syrrunetry  or  trace  properties.  However,  by 
combining  eqs.  (64),  (74),  and  (78),  its  relationship  to  the  F  matrix  can  be  established  as: 

(79) 

Combining  this  result  with  eq.  (35)  we  also  note  that: 

Amv  =  -r^x/M/ve^e^  (80) 
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Since  the  density  p  is  constant,  so  is  the  density  vector  D  .  Consequently,  inserting  eq. 
(78)  into  eq.  (66)  yields  the  following  form  for  the  gravity  gradient  tensor: 

(81) 

Tbe  derivatives  of  the  Q  matrix,  obtained  by  direct  computation  and  subsequent  simplification, 
are  listed  in  Table  S.  When  eq.  (81)  is  combined  with  eqs.  (35)  and  (76),  the  following  basic 
result  is  obtained: 

Gp  Afiv  (82) 

Since  the  gravity  gradient  is  known  to  be  a  tensor  by  virtue  of  the  fact  that  it  is  the  gradient 

of  a  vector,  and  since  the  factor  Gp  is  just  a  constant,  eq.  (82)  implies  that  the  A  transformation 
matrix  encountered  in  the  magnetic  case  is  also  a  tensor  as  is  its  gradient.  Furthermore, 
combining  eqs.  (67a) ,  (68),  and  (82)  leads  to  eq.  (37b). 

A  variety  of  additional  relationships  among  the  various  components  of  the  gradieni  of  die  Q 
matrix  can  be  obtained  by  examining  the  symmetry  properties  of  the  A  tensor  in  the  context  of 
eq.  (35)  and  also  by  direct  examination  of  Table  5.  Several  of  these  relationships  are  listed  in 
Table  6. 

3.4  Gravity  Inversion 

Although  there  are  some  differences,  the  gravitational  inversion  procedure  follows  the  same 
basic  line  of  thought  as  was  outlined  for  the  magnetic  case.  We  first  note  that  multiplication  of 


60 


61 


62 


eq.  (82)  by  ,  the  elements  of  the  inverse  of  A,  and  subsequent  contraction  (i.e.,  setting  two 
indices  equal  to  each  other  causing  a  summation)  of  the  indices  immediately  yields  the  inverse 

relationship: 

Gp  =  (83) 

This  equation  permits  the  determination  of  the  density  p  of  a  single  prism  given  the  geometry 
(i.e.,  the  location  and  dimensions)  of  that  prism  through  the  inverse  matrix  A  and  given 
measurements  of  the  gravity  gradient  field  that  surrounds  the  prism.  However,  we  really  want  to 

be  able  to  determine  not  only  the  density,  but  also  the  geometry  (size  and  location)  of  the  prism. 

This  more  general  problem  is  intrinsically  non-unique  if  all  that  we  have  to  work  with  is  eq. 

(83).  Fortunately,  there  is  more  that  can  be  said. 

Begin  by  denoting  the  inverse  of  the  3  x  3  matrix  Q  as  H  .  The  elements  of  these  two 
matrices  satisfy  the  relation: 

S2|I  Qlv  —  5(iv  (8^) 

wdiere  is  the  Kronecker  delta.  The  elements  of  this  inverse  matrix  are  defined  by  the  relation: 


^  cofactotinvii) 

~  dctCO) 


Then,  quite  formally  the  inverse  of  eq.  (78)  is: 


(85) 


Combining  this  result  with  eq.  (81)  then  gives  the  basic  result: 


(86) 
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9|iv  - 


(87) 


CcHisequently,  a  decoupling  of  the  geometric  parameters  (i.e.,  ,  y,, ,  Zg ,  X, ,  Xy  ,  and  \ )  from 

the  geq>hysical  parameter  (i.e.,  the  density  p)  has  been  accomplished.  The  intention,  thoi,  is  to 
solve  die  nonlinear  problem  indicated  by  eq.  (87)  for  the  geometric  parameters  and  subsequently 
solve  the  linear  problem  indicated  by  eq.  (83)  for  the  geophysical  parameter.  Due  to  the 
symmetry  and  trace  properties  of  the  gravity  gradient  tensor  which  leave  five  independent 
conqionents,  there  are  more  unknown  geometric  parameters,  a  total  of  six,  than  there  are 
equations  to  solve  for  them.  Consequently,  either  one  uses  gravity  measurements  obtained  at 
several  locations,  or  one  may  make  some  restrictive  assumptions  about  either  the  dimoisions  or 
location  of  the  prism.  In  either  case,  the  solution  can  be  obtained  by  stochastic  or  generalized 
inversion  techniques  that  attempt  to  minimize  a  chi>square  (x^)  function  constructed  from  the 
differences  between  the  gravitational  gradients  computed  from  eq.  (87)  and  the  measured 
(observed)  gravity  gradients.  Note,  however,  that  it  is  not  necessary  to  measure  the  gravity 
gradioits  directly  since  diey  may  be  computed  from  a  rectangular  harmonic  potential  function 
doived  from  measurements  of  the  Z-component  of  the  gravity  acceleration.  Also,  advantage 
can  be  taken  of  the  fact  that  the  X-component  and  Y-component  of  the  gravity  acceleration 
vector  can  be  computed  in  terms  of  the  Z-component.  The  general  procedure  is  the  same  as  that 
outlined  for  the  magnetic  survey  of  the  Juan  de  Fuca  region. 

It  is  useful  to  examine  eq.  (78)  in  more  detail.  It  represents  the  following  three  equations: 

gi  =  Gp(D||  11|2  +  Qij)  (88a) 


64 


gj  ss  Gp(fi21  +  Q22  +  ^^23) 


(88b) 


g3  *  GpC^si  +  Q32  +  Q33) 


(88c) 


Taking  the  ratio  of  eq.  (88a)  with  eq.  (88c)  and  taking  the  ratio  of  eq.  (88b)  with  eq.  (88c)  then 
yields: 


gi 


(Q|i  +  0|2  O13) 

(G31  +  032  +  O33) 


(89a) 


(O21  O22  Q23) 

^  “  (03,+  0j2+  033) 


(89b) 


These  equations  show  explicitly  the  dependence  of  the  g,  and  gj  components  of  the  gravitational 
acceleration  vector  on  the  third  component,  gj ,  which  is  the  component  that  is  most  frequently 
measured.  From  eq.  (88c)  we  also  find  that: 

op  =  w 

This  equation  can  be  used  to  determine  the  density  of  the  prism  once  eq.  (87)  has  been  used  to 
determine  the  geometric  parameters  and  hence  £2 . 


3.5  Inverse  Gravity  Modeling  with  Multiple  Prisms 

Generalizing  to  multiple  prisms,  we  now  consider  the  Earth's  crust  beneath  a  gravitationally 
surveyed  area  to  be  foliated  with  a  total  of  N  prisms,  each  of  which  has  its  own  uniform 
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density,  p, .  The  gravitational  acceleration  vector  and  the  gravitational  gradient  tensor  at  the  k'th 
observation  point,  r,^ ,  are  then  composite  sums  of  the  gravity  fields  generated  by  the  individual 
prisms.  Therefore: 


gHk  =  X  Dx.  (91a) 

B-l  '' 

N  .  . 

»Mvk  =  Z  «//v  Dx»  (91b) 

D*l  ^ 

where,  as  in  the  magnetic  case,  the  subscript,  k,  identifies  a  particular  point  of  observation, 
such  that  g^-g^(rj,  fn/]  s  Q/(rk.ri;),  and  [q/zv]  *  «//v(rky.). 

Dj,  is  the  X  component  of  the  density  vector  in  the  n'th  prism,  which  is  centered  about  the  point 
r* 

In  more  succinct  fashion,  eqs.  (91a)  and  (91b)  may  be  written  in  terms  of  matrices  as: 

(92) 


B|1V  —  S|1^vB)l  (93) 

where,  assuming  there  exists  a  total  of  K  observations  for  each  component  of  the  gravity  field 
and  the  gravity  gradient  field,  the  matrices  in  the  above  equations  take  the  following  explicit 
forms: 
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Du 

Du 


(98) 


I,  Din  j 


A  fonnal  inversion  of  eq.  (92)  yields  the  density  matrix  ®  ^  in  terms  of  the  gravity 
acceleration  matrix  6^ : 

Sx  =  Svxe’*''®^  (99) 

where  0**'"  is  the  transpose  of  the  matrix  e»‘''  with  respect  to  the  Latin  indices  n  and  k,  and 
where  the  matrix  Avx  is  the  inverse  of  the  N  x  N  matrix  Ayx,  which  is  defined  as  follows: 

Avi  =  e'‘v0MX  (100) 

Inserting  eq.  (99)  into  eq.  (93)  yields  the  desired  sq>aration  of  the  geometric  and  geophysical 
parameters: 

=  S|i\Apx0  (101) 

The  right-hand  side  of  this  equation  consists  of  27  terms  which  result  from  the  three  Einstein 

summations  over  the  dummy  indices  a,  p,  and  X  .  Choosing  as  an  example,  just  one  of  these 

terms  which  corresponds  to  the  index  specification:  ^  =1,  v  =  2,  X  =  3,  a  =  1,  and  p  =  2,  we 
obtain  the  typical  term:  H,  ’j  A23  0  ,  where  is  a  K  x  N  matrix,  A23  is  an  N  x  N  matrix,  e 

is  an  N  X  K  matrix,  and  61  is  a  K  x  1  matrix. 
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Equation  (101)  is  a  complicated  nonlinear  system  of  equations  involving  6N+1  unknown 
geometric  parameters;  the  number  of  prisms  N,  being  the  one  extra  unknown.  Had  we  not 
already  pre-specified  the  orientation  of  the  prisms  as  being  aligned  parallel  to  the  coordinate 
axes,  there  would  be  an  additional  three  unknown  Eiiler  angles  per  prism  to  be  evaluated, 
yielding  9N+1  unknowns.  Problems  such  as  these  can  be  solved  by  stochastic  inversion  or 
generalized  inversion  techniques.  Such  solutions,  having  been  optimized,  are  therefore  unique 
in  the  least-squares  sense  as  being  the  “best  estimate"  of  the  unknown  parameters.  Once 
determined,  the  geophysical  parameters  follow  directly  from  eq.  (99). 

3.6  The  Computational  Algorithm  GRVREP 

Subroutine  GRVREP  (GRaVity  field  due  to  a  REctangular  Prism)  computes  the  gravity 
potential  U(r)  from  eq.  (74),  the  gravity  acceleration  vector  g(r)  from  eq.  (78),  and  the  nine 
elements  of  the  gravity  gradient  tensor  9(r)  from  eq.  (81)  or  equivalently,  eq.  (82).  It  is 
assumed  that  the  user  of  this  routine  has  written  a  driver  program  for  the  GRVREP  subroutine 
that  defines  the  prism  parameters  and  passes  them  to  GRVREP  through  a  common  block  called 
/GRVBLK/.  It  is  further  assumed  that  the  location  and  orientation  of  the  prism  are  related  to  the 
origin  in  the  same  manner  as  for  Subroutine  MAGREP,  as  discussed  in  Section  2.7.  The 
FORTRAN  code  for  this  subroutine  is  internally  documented  and  is  listed  in  Appendix  F. 

4.  THE  ELECTRIC  FIELD 

Considering  again  only  static  sources,  the  electric  field  can  be  analyzed  using  the  results 
already  obtained  for  the  magnetic  and  gravity  fields.  There  are  two  separate  situations  that  can 
be  considered.  The  first  assumes  that  the  source  of  the  field  consists  of  a  charge  density  a(r') 
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that  is  contained  in  some  volume  V.  The  second  assumes  that  the  source  of  the  field  is  a  material 
within  the  volume  V  containing  an  electric  polarization  P(r').  In  either  situation,  it  is  assumed 

that  there  are  no  electric  currents  and  that  the  volume  V  can  be  broken  into  one  or  more 

subvolumes  In  each  of  which  the  charge  density  and  the  polarization  may  be  considered  uniformly 
distributed. 

4.1  Fields  Generated  by  an  Electrically  Charged  Medium 

The  mathematical  analysis  of  a  static  electrically  charged  medium  has  a  one-to-one 
correspondence  to  that  of  the  gravitational  case  previously  described  since,  under  the 
electrostatic  assumptions.  Maxwell's  equations  reduce  to  the  following: 

V*E  =  47ca  (102a) 

VxE  =  0  (102b) 

These  equations  are  identical  in  form  to  eqs.  (55a)  and  (55b).  Therefore,  the  results  regarding 
gravity  fields  in  section  (3)  may  be  taken  over  completely  if  we  replace  g  by  E  and  -Gp  by  a  . 
We  also  chose  to  make  appropriate  notational  replacements  such  that  U  and 

-»  correspond  to  the  electrostatic  potential  (units:  statvolts),  the  electric  field  gradient 
tensor  (units:  statvolts/cm^),  and  the  electric  charge  density  vector  (units:  statCoulombs/cm’), 
respectively.  Using  thir  notation,  the  results  for  a  single  rectangular  prism  containing  a 
imiform  charge  density  are  summarized  along  with  magnetic  and  gravity  results  in  Table  7. 


Table  7.  Summary  of  Field  Equations  for  the  Single  Prism  Geometry 


GRAVITY 

ELECTRIC 

ELECTRIC 

MAGNETIC 

Source: 

Source: 

Source: 

Source: 

uniform  mass  density 

uniform  charge  density 

uniform  poiarizalion 

uniform  magnetization 

u  =  -  roxe® 

r=  -r^xe^Q" 

T'P 

=  -Q^e“Px 

O  =  -  tl^e-Mx 

u  =  - 

r=  -Foxn^VE^ 

T'P 

=  -Q^Ax'‘£»Epp 

0  =  -  n‘^„Ax‘‘£“Bp 

g4  =  £i/D)L 

=  n/Qx 

Epu 

=  A/Px 

Bp  =  Ap^Mx 

<Sjiv  —  Qjj  /vQx 

^pjiv 

=  A//yPx 

^PV  =  Ajj^/yMx 

^S(iv  =  fln^/ynx,**go 

^jiv  =  Qj,^/yQx  F<j 

^ppv 

=  A/,yAx“Ep« 

^pv  =  A^*'/yAx  Bo 

^fiv  “  GpAjiv 

^jiv  "  “"OAjiv 

< 

II 

< 

^|iv  =  <Sv(l 

<»ptiv 

=  <5pvp 

^pv  =  ^vp 

_  -«iiGp  rOV' 

0  r  >  3v' 

^  _  4«o  r  S9V' 

0  r>dv' 

<?P(1  ; 

=  0  r  >  av' 

^"p  =  0  r  >dV' 

Dx  “  6v 

Qx  =nx'‘E, 

Pk 

=  VEpp 

Mx  =  Ax'*  Bp 

Di=  Gpex 

Qx  =  -  aex 

Gp  =  jA'‘''S^v 

a  =  -  f 

p 

=  yp"Px 

M  =  Jm^Mx 
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4.2  Fields  Generated  by  an  Electrically  Polarized  Medium 

Maxwell's  equations  for  this  case  take  the  following  electrostatic  form; 

V*D  =  0  (103a) 

VxEp  =  0  (103b) 

These  equations  are  identical  in  form  to  eqs.  (4a)  and  (4b)  of  the  magnetic  case.  Consequently, 

the  results  of  that  case  may  be  taken  over  completely  if  we  make  the  appropriate  notational 
rq)lacements  M-^P,  ,  and  .  Here,  is  the  electric  potential  (units: 

statvolts),  and  is  the  electric  field  gradient  tensor  (units:  statvolts/cm^)  due  to  the 

polarization  P  .  The  subscript  (p)  has  been  introduced  to  distinguish  the  present  electrically 
polarized  medium  case  from  the  electrically  charged  medium  case  of  the  previous  section. 
Using  this  notation,  the  results  for  a  single  uniformly  polarized  prism  are  also  summarized  in 
Table  7.  These  seem  to  be  the  least  used  of  the  geopotential  relations  derived  in  this  report. 
However,  it  should  be  noted  that  other  non-geophysics  applications  for  them  do  exist. 

5.  CONNECTIONS  AMONG  THE  ELECTRIC,  MAGNETIC,  AND  GRAVITY  FIELDS 
Over  a  century  and  a  half  ago,  Poisson  (1826)  noticed  that  close  a  relationship  exists 
between  the  magnetic  potential  and  the  gravity  potential  when  the  sources  of  those  potentials  are 
uniformly  distributed  over  the  same  volume.  Similar  relationships  exist  between  the  electric 
field  and  the  gravity  field  potentials,  as  well  as  between  the  electric  and  magnetic  potentials. 
These  relationships  can  be  extended  to  include  the  vector  and  gradient  components  of  these 
fields  as  well. 
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Taking  just  the  relationship  between  the  magnetic  and  gravity  fields  as  an  example,  we  fmd 
that  if  eq.  (31)  is  multiplied  by  the  factor  Gp  and  subsequent  note  of  eq.  (82)  is  taken,  one  of  the 
many  possible  forms  of  Poisson's  relation  is  obtained: 


GpBn  = 


which,  after  multiplication  by  B  ,  gives: 


or,  after  inverting  eq.  (104),  we  have: 


Mx  =  (ix**  Bh)  Gp 


(104) 


(105) 


(106) 


Also,  after  inserting  eqs.  (18),  (64),  and  (66)  into  eq.  (104)  and  subsequently  integrating  with 
respect  to  dx** ,  Poisson's  relation  can  be  expressed  in  standard  form  with  respect  to  the  gravity 
and  magnetic  potentials,  as: 


U,xM^ 


(107) 


where  the  arbitrary  constant  of  integration  has  been  set  to  zero.  Furthermore,  by  taking  the 
gradient  of  eq.  (104)  another  form  of  Poisson's  relation  is  found  to  be: 


GpB|i/y  —  Mx 


(108) 
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All  of  the  above  forms  of  Poisson's  relation  for  the  gravity  and  magnetic  fields  are  expressed  in 
terms  of  the  geophysical  sources  parameters  (i.e.,  magnetization  and  density)  of  those  fields. 
However,  if  eqs.  (106)  and  (108)  are  combined,  then  the  source  parameters  may  be  eliminated. 
The  result  is  a  most  interesting  relationship  between  the  fields  themselves: 

(109) 

Thus,  given  a  priori  knowledge  of  the  gravitational  field,  this  equation  may  be  considered  as  a 
coupled  system  of  linear  Hrst-order  differential  equations  for  the  three  components  of  the 
magnetic  field.  This  can  be  seen  more  clearly  ifeq.  (19)is  inserted  on  the  leff-hand  side  of 
eq.  (109),  thereby  yielding: 

Bk/v  =  (110) 

This  equation  embodies  the  essence  of  a  classical  unification  of  fields  since  it  permits  one  of  the 
two  otherwise  unrelated  fields  (i.e.,  gravity  and  magnetism)  to  be  derived  from  the  other  without 
reference  to  any  additional  fields  or  source  parameters.  This  relationship  goes  well  beyond  the 
standard  Poisson  relations  for  these  fields.  From  this  point  of  view,  it  seems  reasonable  to 
suggest  that  eq.  (110)  and  other  equations  similarly  derived  for  the  electric  field  might  be 
considered  as  classical,  static,  weak-fleld-limit  approximations  of  a  true  classical  unified  field 
theory  (i.e.,  an  extension  of  General  Relativity)  and  that  eq.  (110)  ought  to  be  considered  as 
supporting  evidence  for  the  existence  of  such  a  theory. 

In  defense  of  this  notion,  we  mention  as  an  aside  that  there  is  a  connection  between  the  A 
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matrix  and  certain  components  of  the  Riemann  curvature  tensor,  ^  ,  of  General  Relativity. 

The  Riemann  curvature  tensor  is  a  tensor  of  rank  four,  which  characterizes  the  geometry  of  a 
4-dimensional  Riemann  space.  Here,  Latin  indices  range  from  0  to  3,  while  Greek  indices  still 
range  from  1  to  3.  In  a  4-dimensional  Riemann  space,  the  zero'th  coordinate  is  related  to  time, 
which,  when  multiplied  by  the  speed  of  light,  c,  has  units  of  length.  Then,  in  addition  to  the 
three  space  coordinates  x'  =  x  ,  x^  =  y,  and  x*  =  z  ,  we  also  have  x®  =  c  t .  Details  regarding 
the  definition  of  the  Riemann  curvature  tensor  and  its  relationship  with  the  gravitational  field 
equations  of  General  Relativity  and  its  Newtonian  approximation  are  discussed  by  Adler,  Bazin, 
and  Schiffer  (1975).  Here  we  simply  note  that  linearization  of  the  gravitational  field  equations 
of  General  Relativity  in  the  static,  weak-fleld  approximation  yields  the  following  result: 


po  = 


_  -L  Tj/» 

c2  ^  n 


(111) 


which  establishes  the  connection  between  the  Riemann  curvature  tensor  and  the  A  matrix. 
Furthermore,  if  we  denote  the  inverse  of  the  3  x  3  matrix  ^“®po  as  ^“®po  so  that: 


^’'"«o^“"po  =  8’'p  (112) 

we  find  that  eq.  (110)  can  be  put  into  the  form: 

B^/v  (113) 


A  similar  result  can  be  obtained  for  the  static  electric  field.  Using  results  such  as  these,  it  may 
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be  possible  to  extend  these  field  relationships  to  covariant  4-dimensional  form  and  thereby 
obtain  a  system  of  classical  unified  field  equations,  which  in  the  limit  of  weak  electromagnetic 
fidds  reduce  to  the  tensor  form  of  Maxwell's  equations  and  the  equations  of  Goieral  Relativity. 
That  »ich  a  classical  unified  field  theory  should  exist,  regardless  of  what  further  unifications 
may  exist  at  the  quantum  level,  is  discussed  by  Tonnelat  (1966a,  1966b). 

Several  other  Poisson  relationships  exist  between  the  electric  and  magnetic  fields  and 
between  the  electric  and  gravitational  fields.  They  can  be  derived  almost  by  inspection  in  a 
manner  that  is  essentially  the  same  as  that  described  above  between  the  magnetic  and 
gravitational  fields.  For  one  more  example,  take  the  case  of  the  electric  field  E^(r)  due  to  the 
polarization  P(r)  and  their  relationship  to  the  magnetic  induction  B(r)  due  to  the  magnetization 
M(r).  Using  the  correspondences  established  in  Section  4.2,  we  can  say  that: 

EpH=A/Pi  (114) 

This  has  the  following  inverse  relationship: 

Px  =  A/Ep^  (115) 

Taking  the  vector  dot  product  between  the  polarization  vector  and  the  magnetic  induction  vector 
using  eqs.  ^,1 15)  and  (34)  and  subsequently  taking  note  of  the  symmetry  of  the  A  matrix  and 
using  eq.  (40),  we  obtain  the  following  Poisson  relation: 

B,.P'‘  =  MxE^  (116) 
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Taking  the  gradient  of  this  equation  and  noting  the  definition  of  the  magnetic  gradient  tensor 
given  by  eq.  (19)  and  that  the  electric  field  gradient  due  to  polarization  d‘|,(r)  is  similarly  defmed 

in  tenns  of  the  electric  field  due  to  polarization  Ep(r),  we  find,  due  to  the  uniformity  of  the 
magnetization  vector  and  the  polarization  vector,  another  Poisson  relation: 

M»-  (117) 

By  multiplying  through  eq.  (1 17)  by  the  inverse  of  the  magnetic  gradient  tensor,  we  find  that: 

Pv  =  (118) 

Alternatively,  multiplying  eq.  (1 17)  by  the  inverse  of  the  electric  field  gradient  tensor,  we  find 
that: 


Mx  =  ^v^JxPp  (119) 

The  symmetry  property  of  the  two  gradient  tensors  was  also  used  to  derive  eqs.  (118)  and 
(119).  These  and  other  Poisson  relations  are  listed  in  Table  8  and  Table  9.  However,  this 
list  of  relations  is  not  an  exhaustive  one.  Furthermore,  they  are  not  exclusively  applicable  to 
geophysics  problems.  They  can  also  be  applied  to  problems  in  biophysics,  engineering,  etc. 
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Hibie  8.  Poisson  Relations  Among  the  Gravity,  Electric,  and  Magnetic  Fields  1 


Sources:  Sources:  Sources:  Sources: 

mass  density  mass  density  charge  density  charge  density 

&  &  &  & 
polarization  magnetization  polarization  magnetization 


TaUe9. 


Poisson  Relations  Annong  the  Gravity,  Electric,  and  Magnetic  Fields  II 


Sources: 

mass  density 
& 

charge  density 

Sources: 

magnetization 

& 

polarization 

Gp  =  -  f  ?*'' 

Ml  = 

0= 

Pm  = 

C 

II 

1 

o  —  _  ^  F 

SM  0 

P^B^  =  EpiM" 

—  2L  p 

cr|iv  —  “  C?|iv 

PI'^^v  =  ^pvxM*- 
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APPENDIX  A 


DERIVATION  OF  THE  Q  MATRIX  ELEMENTS 


Recalling  eq.  (12),  which  is  an  expression  for  the  magnetic  scalar  potential  <l>(r)  in  terms  of 
the  vector  magnetization  M(r').  we  now  specify  that  this  magnetization  be  uniformly  distributed 
over  the  volume  V  of  a  single  rectangular  prism  and  hence  is  a  constant  vector  M,  so  that: 


«l>(r)  = 


(Al) 


In  Cartesian  coordinates  this  equation  takes  the  form: 


<I»(X,y,Z)  =  Mxlx  +  Myly  +  Milz 


(A2) 


where  M, ,  M^, ,  and  are  the  constant  components  ot  the  magnetization  vector  M,  and  where 
the  following  notation  has  been  used: 


i*{x.y,z) 


iy(x.y,z) 


Iz(x,y,2) 


dx'  dy^dz' 


-f- 


dx'  dy'  dz' 


dx'dy'dz' 


(A3a) 


(A3b) 


(A3c) 


where; 


A-3 


R  *  I r  -  r'  I  =  J^x  -  x'j  +  (y  -  /]  +  -  z'J  (A4) 

The  integrals  I,  ,  I,  ,  and  are  just  cyclic  permutations  of  the  same  mathematical  form. 
Therefore,  it  is  sufficient  to  evaluate  just  one  of  these  integrals,  for  instance  1, ,  and  obtain  the 
others  by  permuting  x,  y  and  z  appropriately  wherever  they  occur.  Thus,  taking  1^  and 
performing  the  indicated  derivative  with  respect  to  x  gives: 

I,(x,y,z)  =  J  ^  dx'd/dz'  (A5) 

The  volume  of  integration,  V  ,  is  taken  to  be  a  rectangular  prism  centered  about  the  point 
ffl  =  (Xfl  ,yo,Zo)  with  its  sides  parallel  to  the  coordinate  axes  and  having  dimensions  \ ,  and 
X^ .  Consequently, 


zy  yv 

I«(x,y,z)  =  J  J  J  dx'd/dz' 

zl  Yl  xl 


where  the  upper  and  lower  limits  of  integration  are  defined  as  follows: 

A.* 

XL  =  Xo  -  — 

Xx 

Xy  =  Xo  +  — 


(A6) 


(A7a) 

(A7b) 
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yt  =  yo  -  -y 


(A7c) 


yu 


(A7d) 


Zl  =  Zo  - 


2 


(A7e) 


zu  =  zo  +  T 


(A7f) 


In  the  event  that  we  were  dealing  with  multiple  prisms,  each  term  in  the  above  set  of  equations 
would  include  the  additional  subscript  n,  which  identifies  a  particular  prism.  Integrating  first 
with  respect  to  x' ,  gives: 


zu  yu 

Iz(z.y,z)  =  J  J  lx/ 

zl  XL 


(A8) 


where  we  have  used  the  shorthand  notation: 


x'  =  xu 
x'  =  Xl 


(A9) 


to  indicate  the  limits  of  evaluation  of  the  x'  coordinate.  Subsequently,  the  same  notation  will 
also  be  used  for  the  limits  of  the  y’  and  z’  integrations. 

Next,  integrating  over  the  y'  coordinate,  the  following  result  is  obtained: 


A-5 


Zu 


I»(x,y,z)  =  -  J  ln[(y  -  yO  +  R]dz'  Ij/.y/ 


(AlO) 


The  first  two  integrals  over  x'  and  y'  are  found  in  standard  mathematics  tables  of  integrals.  The 
last  integral  over  z'  is  not  as  easy  to  evaluate.  In  order  to  perform  this  last  integration,  first 
define; 

u  =  (y  -  /)  +  R  (A1 1) 

Using  eq.  (A4)  for  R,  the  above  equation  can  be  solved  for  (z  -  z')  in  terms  of  the  variable  u: 
z  -  z'  =  ±  S(u)  (A12) 

where; 


S(u)  = 


-  (x  -  x')  -  (y  -  /) 


(A  13) 


Taking  the  derivative  of  eq.  (A  12)  then  gives: 

dz'  =  ±  dS(u)  (A  14) 

Now,  S(u)  is  always  positive  since  it  is  just  equal  to  |  z  -  z'  | .  So,  the  ±  signs  above  must  be 
chosen  to  correspond  with  the  sign  of  (z-z') .  Thus,  the  positive  sign  (+)  is  chosen  if 
z  -z'  s  0  ,  while  the  negative  sign  (-)  is  chosen  if  z  -  z'  <  0  ,  the  sign  changing  at  the  point 
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2  =  2'.  Consequently,  the  integral  is  broken  into  two  parts  as  follows; 


I*(x.y.z)  =  -  j  In[(y-y']  +  Rjdz'  +  ln[^(y  -  y' j  +  Rjdz'j  Ix'.y'  (A15) 


Then,  in  terms  of  the  variable  u,  the  integrals  in  eq.  (A  15)  may  be  written  in  the  following 


manner; 


IfuU)  fUu  I 

In(u)dS(u)  +  ln(u)dS(u)  Ix'y' 
JuL  Ju(z)  ■' 


(A16) 


where,  by  the  notation  and  we  mean  u(z*  =  zj  and  u(z'  =  z^)  respectively.  Integrating 
by  parts  gives; 


Ix(x,y.z)=  -  -  S(u)  ln(u)  +  S(u)  ln(u)  ^  ^du  -  T  -^J^i-dul  Ix/.y/  (A  17) 

JuL  Ju(l)  J 


The  argument  underneath  the  square  root  in  S(u)  is  just  quadratic  in  u.  Therefore,  the  two 
integrals  in  eq.  (A17)  may  be  partially  evaluated  using  Gradshteyn  and  Ryzhik  (2.2671), 
yielding; 


Ix(x,y,z)  =  -  -  S(u)  ln(u)  lllf’  +  S(u)  ln(u)  1"^",  +  S(u)  IHf’  -  S(u)  -  (x-x')^  P* 

^  ul 

,  r**(2)  ruu  1  I 

"  L  ^  ^ L 


Finally,  using  Gradshteyn  and  Rhyzik  (1980)  sections  2.261  and  2.266  it  is  found  that; 
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I,(x.y,z)=  -|[l-In(u)]S(u)li*^  -  [1 -ln(u)]S(u)  -  (x-x')  tan"' 

+  (x-x')  Un-'f  ->->>]  IS,  -  (,  -  y-)  totsu)  +  X  -  (y  -  y-)]  C 


+  (y  -  yO  ln[S(u)  +  u  -  (y  -  /)]  I x/ ,  y/ 


Recalling  the  definition  of  u  given  in  eq.  (AlO),  the  above  equation  becomes: 


lK(x,y,z)  =  -  (Iz  -  z'l  { I  -  ln[R  +  (y  -  /)]}  -  Iz  -  z'l  {I  -  In[R  +  (y  -  /)]}  IJ" 

X..  ..A  x-.-if -(»-»')’ -(y-yV-{y-y')Rl  u  .  x..  ..a _ if -(*-»V-(y-yV-(y-y')R  1  nu 

-(x-xOtan'l^ - - J - (,-,/) - J 


-(y-/)  ln[R+ |z-z'|]liL  +  (y-yO  ln[R  +  |z-z'|]lz“}  Ix'.y' 


This  result  simplifies  somewhat  when  it  is  noted  that  the  following  identity  holds: 


tan"'(-Q)  =  -  tan"'(Q) 


(A19) 


(A20) 


(A21) 


Then,  when  the  factor  |  z  -  z'  |  is  replaced  by  ±(z  -  z')  depending  on  the  limits  of  the  z' 
int^ration,  it  is  found  that: 

i.(x.y,x)-  -((z-x'){i-  u.[R+(y-y')]iis  (x  - .')  m-'[  ]  la 


-(y-yO  InCR  +  (z-z')]Ih  +  (y  -  y')  ln[R  -  (z  -  z')  I^]}  Ix'.y' 


(A22) 
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Further  note  that  any  term  that  does  not  involve  all  three  coordinates  x' ,  y' ,  and  z*  will,  upon 
evaluation  at  the  limits  of  integration,  be  equal  to  zero.  Consequently,  the  first  term  in  eq.  (A22) 
IS  zero.  That  is: 


(z— z^)  Ix'.y^.z^  =0 


(A23) 


Using  the  identity: 


InfR  —  (z  —  z03  Ix^.y^.z^  —  —  ln[R  +  (z  — z03  lx^,y^,z^ 


(A24) 


the  validity  of  which  depends  on  the  fact  that  these  logarithms  are  being  evaluated  at  fixed 
limits,  I,(x,y,z)  reduces  to: 

I*(x,y,z)  =  -  j(x-x')  -  (z  -  t!)  ln[R  +  (y  -  /)] 


-(y-/)ln[R  +  (z-z')]}  Ix'.y'.z' 


(A25) 


Next,  applying  the  identity: 


[tan-'(Q)  +  tan“'(Q-')]  L'.y'.z'  =  £  f  L'.y'.z'  =  0 


(A26) 


where: 


+1  Q  >  0 
-I  Q  <  0 


(A27) 
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and  subsequently  performing  some  straightforward  manipulations  to  simplify  the  result,  we  find: 


(x-»')(z-z')  1  \  , 


(y-y')U-z')  1 

(x-»')R  J 


(A28) 


Consequently,  the  final  result  is  obtained: 


Ix(x,y,z)  =  - 1+  (X  -  x')  tan-'  [  -  ]  -  (y-/)  ln[R  +  (z-z')]  -  (z-z')  ln[R  +  (y-/)]J  Ix'.y'.z' 


(A29a) 


The  integrals  Iy(x,y,2)  and  I/x,y,2)  can  be  obtained  from  the  above  expression  for  I^(x,y,2)  by 
cyclic  permutations  of  the  coordinate  differences,  i.e.: 

(x-xO  (y-/)  (2-2O  -»(x-xO 

yielding: 


Iy(x,y,z)  =  -  |-(x-xO  ln[R  +  (z-zO]  +  (y-y')fan  '[  - zO  In [R  +  (x - xO]}  lx^y^z' 

(A29b) 


and 


Iz(x,y,z)  =  -  |-(x-x')  ln[R  +  (y-/)]  -  (y-y')  In[R  +  (x-xO]  +  (z- z^  Ix'.y'.z^ 


(A29c) 
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The  magnetic  scalar  potential  is  therefore  determined  and  is: 


0(x,y,z)  =  M.I>(x,y,z)  +  Myly(x,y.z)  +  MJ,{x.y,z)  (A30) 

The  elements  of  the  Q  matrix,  listed  in  Table  1,  are  defined  to  be  the  nine  terms  contained  in  the 
three  functions  I,(x,y,z),  Iy(x,y,z)  and  ljj(x,y,z).  Defining  the  matrix  elements  in  this  way  and 
given  the  definition  of  the  vector  e,  eq.  (A30)  can  be  put  into  the  form  of  eq.  (30)  of  the  main 
text 


All 
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DERIVATION  OF  THE  A  MATRIX  ELEMENTS 
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In  terms  of  the  magnetic  scalar  potential  <I>(x.y,z)  the  magnetic  field  vector  B(x,y,z), 
generated  by  a  single  prism,  is: 

B(x.y.z)  =  -  V<I)(x.y.z)  (Bl) 

Using  the  form  for  the  scalar  potential  given  in  eq.  (A2)  with  given  in  the  form  of  eq.  (A  10), 
while  ly  and  I,  are  obtained  from  eq.  (A  10)  by  cyclic  permutation  of  the  coordinate  differences 
as  previously  described  in  Appendix  A,  the  magnetic  induction  B  can  be  cast  into  the  following 
form: 


B(X,y,Z)  =  Mk.#X  +  My,^y  +  MiJz 


where  the  vectors  ,./y,and,/j  are  defined  as  follows: 


*^*(x.y.z)  =  f  Vln[(y-/)  +  Rjdz'  Ix'.y/ 

Jyix,y,z)  =  f  “  Vln[(z-z')  +  R]dx'  ly^z/ 

Jz(x,y,z)  =  f  "  V  lii[(x  -  x')  +  R]d/  Ix'.z/ 

JyL 


(B2) 


(B3a) 

(B3b) 

(B3c) 


Since  the  vectors  and  are  all  of  the  same  mathematical  form,  it  is  sufficient  to 

evaluate  just  one  of  them,  for  instance  ,  and  then  to  obtain  expressions  for  the  other  two 
vectors,  via  cyclic  permutation  of  the  coordinate  differences  as  was  done  in  Appendix  A. 


Additionally,  since  we  are  now  dealing  with  vectors,  a  simultaneous  cyclic  permutation  of  the 
unit  vectors  i ,  j  and  k  is  also  necessary.  That  is: 


(x  -  x')  (y  -  yO  (z  -  z')  — >  (x  -  x') 

i  — ^  j  — >  k  — >  i 


Now,  the  gradient  of  the  natural  logarithm  in  eq.  (B3a)  is  easily  computed,  yielding: 


=  r  ]‘  -  (7)i 


f  _ (y-yO(z-z^)  ]  I  /  I 

I  [(x-x')'  ♦U-xO']  [(x-x')»+(z-z')2]R  J  “j  *“  '  X  .  y' 


(B4) 


Noticing  that  the  first  term  of  the  integral  will  contribute  nothing  to  the  final  result,  when 
evaluated  at  the  limits,  since  it  is  independent  of  the  y'  coordinate  and  further  noticing  that  the 
same  is  true  of  the  fourth  term,  (x,y,z)  simplifies  to: 


^x(x,y,z) 


_ <y-y,') _ 

[(x-x')^  +  (y-yO']R 


r 

»  Zi 


(x-x')(z-z') 

[(x-x02  +  (z-x')']R 


dz'  k 


lx'  y.z' 


(B5) 


These  final  three  integrals  may  be  evaluated  using  Gradshteyn  and  Ryzhik  (2.124,  2.261  and 
2.284),  yielding: 


(B6) 
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The  last  term  can  be  simplified  by  multiplying  the  numerator  and  denominator  of  the  logarithm 
argument  by  the  factor  (y  -  y')  R  •  Then,  the  logarithm  can  be  broken  into  the  sum  of  two 
terms,  one  of  which  will  be  independent  of  y'  and  therefore  will  contribute  nothing  when 
evaluated  at  the  limits.  In  this  way,  the  following  identity  (and  cyclic  permutations  thereof)  is 
obtained; 

i’  *“[  (y-'y^  +  fi]  +  (y-yO]  Ix'.y'.z'  (B7) 

Consequently,  finally  reduces  to: 

»^«(x,y.z)  =  -  ln[R  +  (z  -  z')]i  -  ln[R  +  (y  - /)]  k|  lx^y^2/  (B8a) 

Cyclic  permutations  then  yield  similar  e^cpressions  for  and  t^(x,y,z): 

,^y(x,y,z)  =  {-  ln[P+(z-zO]i  +  tan'* [ ^ ]  j  -  ln[R  +  (x  -  x')]  kj  Ix'.y'.z'  (B8b) 
and 

.^z(x,y,z)  =  j-  ln[R  +  (y- y')]i  -  ln[R  +  (x  -  x')] ]  +  tan-'[ ■]  k]  Ix'.y'.z'  (B8c) 

The  nine  terms  in  the  three  eqs.  (B8a),  (B8b),  and  (B8c),  each  evaluated  independently  at 
the  prism  limits,  form  the  nine  elements  of  the  A  matrix  listed  in  Table  2.  These  matrix 


elements,  so  defined,  allow  eq.  (B2)  to  be  written  in  the  form  of  eq.  (34)  of  the  main  text. 
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DERIVATION  OF  THE  F  MATRIX  ELEMENTS 


If  the  mass  density,  p  ,  is  specified  to  be  uniform  over  the  volume  of  integration  V  and 


consequently  is  independent  of  r',  then  eq.  (58)  for  the  gravitational  potential  U(r)  may  be 
written  as: 

U(r)  =  -Op  1^,  (Cl) 

Taking  the  volume  of  integration  V,  as  in  Appendix  A,  to  be  a  rectangular  prism,  centered  about 
the  point  1*0  “  «  yo  >  Zq  )  having  dimensions  X, ,  Xy ,  and  X^ ,  the  above  equation  may  be 

written  as; 


U(x.y,z)  =  -Gpl(x,y,z) 


(C2) 


where: 


i(x.y,z) 


du'dy'ib' 

R 


(C3) 


and  where  the  limits  of  integration  are  the  same  as  those  defined  in  Appendix  A.  The  first 
integral  over  x*  is  found  in  most  standard  mathematics  tables  of  integrals  and  yields  the 
following  result; 


nyv  I 

ln[(x  -  x')  +  Rjdy'dz'  lx' 


(C4) 


The  integral  over  y'  is  of  the  same  mathematical  form  as  that  of  eq.  (A  10),  which,  after  some 
effort,  was  evaluated,  yielding  Iy(x,y,2)  in  Appendix  A.  Consequently,  we  have: 
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K*.y.*) 


(x-xoj*"  ln[R  +  (y-/)]dz'  +  (y-  lii[R  +  (x-x')]dz' 


(C5) 


The  first  two  integrals  over  z*  in  the  above  expression  are  of  the  same  form  as  the  one  just 
peaformed  over  y' .  So  we  have  immediately: 


I(x,y.2)  =  [(x  - xO  [(X - xO  tan-*[  -  (y-/)  ln[R  +  (z-zO]  -  (z-z^)  ln[R  +  (y-y')])  I*/ 


+  (y-/)(-(x-x')In[R  +  (z-z')]  -  (z-z'lInCR  +  Cx-x')]  +  (y/) I*' 

-  Cz-z')tan-'[^2^^^]d2'j  Ix/.y/  (C6a) 


The  final  integral  in  eq.  (C6a)  is  rather  formidable.  However,  there  is  a  little  trick  that  can 
be  used.  Note  that  the  order  of  integration  of  eq.  (C3)  is  arbitrary  because  the  limits  of 
integration  are  fixed.  So  instead  of  integrating  first  over  x'  and  then  over  y' ,  as  was  just  done, 
we  could  have  integrated  first  over  y'  and  then  over  z’,  or  alternatively,  we  could  have  integrated 
first  over  x'  and  then  over  z*.  The  results  of  these  alternative  orders  of  integrations  is  obtained 
from  eq.  (C6a)  by  cyclic  permutation  of  the  coordinate  differences,  yielding: 
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K«.y..)  -  |-  j;“  *  (y-^)jo,-/).an-'[ii^^] 


I 


-(z-2')ln[R  +  (x-x')]  -  (x-xOln[R  +  (2-z')]}  l»'  +  (z-zO{-(y-/)  ln[R  +  (x-x')] 
-(x-x')ln[R  +  (y-/)]  +  (z-z') |,/[  ly/,2'  (C6b) 

and 

I(x,y,z)  =  ((x-xOj-(z-z')ln[R  +  (y-yO]  -  (y-yO  ln[R  +  (z-z')]  +  (x - x^  tan"' [  ]}  ly' 

-  I’"  (y-/)u.-'[ii^]dy'  *  (z-/)((2-.')a«-'[^=g2gl] 

-  (X  -  x')  In[R  +  (y -  y')]  -  (y  -  yO  In [R  +  (x -  x')]}  ly/ }  I  x' ,  z'  (C6c) 

The  three  forms  of  I(x,y,z)  given  in  eqs.  (C6a),  (C6b),  and  (C6c)  must  be  equivalent.  Therefore, 
by  comparing  these  three  equations,  it  may  be  inferred  that ; 


J*“(x-x0tan-'[^^g3gl]dx'  ly/.2/  =-((x-x')({x-x')tan-'[iig3^] 


-(y-/)In[R  +  (z-z')]  -  (z-z')ln[R +  (y-y')]}}  Ix'.y'.z/  (C7) 
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Cyclic  permutations  of  this  equation  are  also  valid.  Thus,  when  eq.  (Cl)  and  its  permutations 
with  respect  to  the  coordinate  differences  are  inserted  into  equations  (C6a),  (C6b),  and  (C6c),  it 
is  found  that  the  gravitational  potential  of  a  rectangular  prism  takes  the  following  final  form: 


U(x,y.z)  =  - Gp |(x - x') |(x - x')  tan"' [ ■■  ]  -  (y-y')  ln[R  +  (z-z')]  -  (z-z')  ln[R  +  (y-/)]J 
+  (y/lj-tx-xO  ln[R  +  (z-zO]  +  (y-y')tan-'[^ij||i^]  -  (z-z^  ln[R  +  (x-x')]} 

+  (z-z')(-(x-x')ln[R  +  (y-y')]  -  (y-y')  ln[R  +  (x-x')]  +  (z-z')  Ix'.y/.z' 

(C8) 

The  elements  of  the  F  matrix  listed  in  Table  4  are  defined  to  be  the  nine  geometric  terms 
contained  in  eq.  (C8).  Defining  the  matrix  elements  in  this  way  allows  eq.  (C8)  to  be  written  in 
the  compact  form  of  eq.  (74)  of  the  main  text. 
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APPENDIX  D 


DERIVATION  OF  THE  GRAVITY  ACCELERATION  VECTOR  COMPONENTS 


The  gravitational  potential  generated  by  a  rectangular  prism  with  a  uniform  mass  density  p  is 


of  the  form: 


U(x,y,z)  =  -  Gpl(x,y,z) 


where,  as  in  Appendix  C,  we  have: 


l(x,y.z) 


dx'Vdz' 

R 


(Dl) 


(D2) 


and  where  the  limits  of  integration  are  the  same  as  those  defined  in  Appendix  A. 

The  gravitational  acceleration  vector,  g,  is  the  negative  gradient  of  the  potential,  so  that: 

g(x,y,z)  =  -  VU(x,y,z)  (D3) 

So  if  the  x'  integration  in  eq.  (D2)  is  performed  first,  then  the  gravitational  acceleration  takes  the 
following  form: 

r*t/  fXu 

g(x,y,z)  =  -  Gp  V  ln[R  +  (x-x')]d/dz'  1,'  (D4a) 


On  the  other  hand,  since  the  limits  of  integration  are  fixed,  the  y'  integration  could  have  been 
performed  first,  giving: 


8(x.y.2) 


Gp  f  f  V  ln[R  +  (y-y')]dx'dz'  ly 

J  Zl  ^  H 


(D4b) 


Anodier  possibility  is  to  perform  the  z*  integration  first,  which  then  yields: 


*(Xty.*)  =  -  Gp  f  f  V  ln[R  +  (z-z')]dx'd/  I*/  (D4c) 

These  last  three  equations  are  just  cyclic  permutations  with  respect  to  the  coordinate  differences 
of  the  same  basic  mathematical  form.  Furthermore,  note  that  the  y'  integration  of  eq.  (IMa)  is 
the  same  as  the  one  previously  encountered  in  eq.  (B3c),  while  the  z'  integration  in  eq.  (D4b)  is 
the  same  as  that  already  encountered  in  eq.  (B3a).  Likewise,  the  x'  integration  in  eq.  (D4c)  is 
the  same  as  dtat  previously  encountered  in  eq.  (B3b).  The  results  of  these  three  integrations  are 
given  in  eqs.  (B8a),  (B8b),  and  (B8c)  of  Appendix  B,  so  that  eqs.  (D4a),  (D4b),  and  (D4c)  can 
immediately  be  written,  respectively,  as: 

g(x,y,z)  =  -  Gp  {-  ln[R  +  (y-/)]i  -  ln[R  +  (x-x')]j  +  tan'* kid/  1,/,,^  (D5a) 
g(x,y,z)  =  -  Gp  I  “  {+*»“■'[ -  In[R  +  (z-z')]j  -  ln[R  +  (y-y')]k|dx'  \y> (D5b) 
g(x,y,z)  =  -  Gp  |- ln[R  +  (z-z')]!  +  -  ln[R  +  (x - x')] k} dy'  I,/.*/  (D5c) 

These  three  expressions  for  the  gravitational  acceleration  vector  must  be  equal.  Therefore,  by 
conqtaring  like  components  of  these  expressions,  we  find  the  following  set  of  identities: 

J"  ^  in[R  +  (2_z/)jdy/  In[R  +  (y-y')]dz^ 

(D6) 
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Cyclic  permutations  of  these  relations  with  respect  to  the  coordinate  differences  are  also  valid. 
The  reason  for  going  through  the  above  exercise  is  that  the  integrations  involving  the 
arc-tangents  are  very  difficult  to  evaluate,  unless  they  are  cast  into  some  other  equivalent  forms, 
which  has  now  been  done. 

Picking  any  one  of  the  three  forms  of  g(x,y,z),  for  instance  that  of  eq.  (D5c),  and  using  one 
of  the  identities  in  eq.  (D6)  for  the  arc-tangent,  we  may  write  the  gravitational  acceleration 
vector  in  the  following  form: 


g(x,y,2)  =  Gpjj^  ln[R  +  (z-zO]dy'  Ix'.z' i  +  J  ln[R  +  (x-x')]dz'  Ix'.y' j 


+  f  ln[R  +  {x-x')]dy'  Ix'.z' k 
•'yt 


(D7) 


All  of  the  integrals  in  this  equation  are  of  the  same  basic  mathematical  form  as  eq.  (A  10)  of 
Appendix  A.  Consequently,  when  evaluated,  the  above  integrations  will  yield  forms  of  the  type 
given  in  eqs.  (A30a),  (A30b),  and  (A30c),  the  particular  form  depending  on  which  permutation 
of  the  coordinate  differences  is  involved.  In  this  way,  the  following  result  is  obtained: 


g(x,y,z)  =  gxi  +  gyj  +  gzk 


(D8) 


where: 

gx(x,y,z)  =  Gp {+  (X - xO tan"' [ -  (y-/)ln[R  +  (z-z')]  -  (z-z')ln[R  +  (y-/)]| lx' .y' ,z' 


gy(x,y,z)  =  Gp(-(x-x')ln[R  +  (z-z')]  +  (y - /) tan"' [ ]  -  (z-z')ln[R  +  (x-x')])  Ix'.y' ,z' 
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r’ 


gi(x.y.z)  =  Gpj-(x-x')ln[R  +  (y-y')]  -  (y - yO In [ R  +  (x - xO ]  +  (z - z') tan~‘ [ H ^ 


Using  the  elements  of  the  fl  matrix  derived  in  Appendix  A  and  listed  in  Table  1,  the  vector 
components  of  the  gravitational  acceleration  given  in  eq.  (D8)  can  be  cast  into  the  form  of  eq. 
(78)  of  the  main  text. 
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THE  MAGREP  FORTRAN  ALGORITHM 
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C************************^'*************'k>‘******************************* 

c 

c 

C  SUBROUTINE  MAGREP  (MAGNETIC  FIELD  DUE  TO  A  RECTANGULAR  PRISM) 

C 

c 

c*********************************************************************** 

c 

c 

c  PROGRAMMED  BY;  JOHN  M.  QUINN  AND  DONALD  L.  SHIEL  8/8/86 

C  NAVAL  OCEANOGRAPHIC  OFFICE 

C  STENNIS  SPACE  CENTER,  MS  39522-5001 

C 
C 

C* ********************************************************************** 

C 

C 

C  PURPOSE:  THIS  ROUTINE  COMPUTES  THE  MAGNETIC  FIELD  COMPONENTS  DUE 

C  TO  A  RECTANGULAR  PRISM  LOCATED  AT  (XB,YB,ZB)  RELATIVE  TO 

C  SOME  ORIGIN  LOCATED  AT  THE  OCEAN  SURFACE  AT  THE  LOWER 

C  LEFT-HAND  CORNER  OF  THE  SURVEY  AREA.  THE  PRISM  HAS 

C  DIMENSIONS  (LAMX,LAMY,LAMZ) .  ITS  ORIENTATION  IS 

C  DESCRIBED  BY  EULER  ANGLES  (ALPHA, BETA, GAMA) 

C  CORRESPONDING  TO  YAW,  PITCH,  AND  ROLL  ACCORDING  TO  THE 

C  3-2-1  CONVENTION,  RELATIVE  TO  THE  USUAL  GEODETIC 

C  COORDINATES  FOR  WHICH  X=NORTH,  Y=EAST,  AND  Z=DOWN.  THE 

C  PRISM  HAS  UNIFORM  MAGNETIZATION  (MX,MY,MZ)  WITH  RESPECT 

C  TO  THE  PRISM-FIXED  COORDINATES. 

C 

C 

c******* **************************************************************** 


C 

c 

c 

c 

C 

c 

C 

c*** 

C 

C 

c 

C 

C 

c 

c 

c 

c 

C 

C 

C 

c 

C 

C 

C 

C 

c 

C 

c 

C 


REFERENCE:  A  UNIFIED  APPROACH  TO  GEOPOTENTIAL  FIELD  MODELING; 

BY:  JOHN  M.  QUINN  AND  DONALD  L.  SHIEL;  U.  S.  NAVAL 
OCEANOGRAPHIC  OFFICE  TECHNICAL  REPORT  No.  308  (1993) 


******************************************************************** 


PARAMETER  DESCRIPTIONS: 

X  -  INERTIAL  X  (NORTH)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

Y  -  INERTIAL  Y  (EAST)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

Z  -  INERTIAL  Z  (DOWN)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

XB  -  INERTIAL  X  (NORTH)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

YB  -  INERTIAL  Y  (EAST)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

ZB  -  INERTIAL  Z  (DOWN)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

XP  -  PRISM  FIXED  X-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

YP  -  PRISM  FIXED  Y-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

ZP  -  PRISM  FIXED  Z-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

LAMX  -  DIMENSION  OF  PRISM  ALONG  X-AXIS  (KM) 

LAMY  -  DIMENSION  OF  PRISM  ALONG  Y-AXIS  (KM) 

LAMZ  -  DIMENSION  OF  PRISM  ALONG  Z-AXIS  (KM) 

ALPHA  -  YAW  ROTATION  ANGLE  ABOUT  Z-XIS  OF  PRISM  (DEG.) 

BETA  -  PITCH  ROTATION  ANGLE  ABOUT  Y-AXIS  OF  PRISM  (DEG.) 

GAMA  -  ROLL  ROTATION  ANGLE  ABOUT  X-AXIS  OF  PRISM  (DEG.) 

MX  -  X-COMPONENT  OF  MAGNETIZATION  OF  PRISM  (NT) 

MY  -  Y-COMPONENT  OF  MAGNETIZATION  OF  PRISM  (NT) 
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MZ  -  Z-COMPONENT  OF  MAGNETIZATION  OF  PRISM  (NT) 

ROT  -  ROTATION  MATRIX  FROM  INERTIAL  TO  PRISM  FIXED  COORDS. 
ROTI  -  INVERSE  ROT.  MATRIX  —  PRISM-FIXED  TO  INERTIAL  COORDS. 
V  -  GEOMAGNETIC  POTENTIAL  IN  INERTIAL  COORDINATES  (NT-KM) 

BX  -  INERTIAL  X-COMPONENT  OF  OBSERVED  MAGNETIC  FIELD  (NT) 

BY  -  INERTIAL  Y-Cc  >ONENT  OF  OBSERVED  MAGNETIC  FIELD  (NT) 

BZ  -  INERTIAL  Z-COMPONENT  OF  OBSERVED  MAGNETIC  FIELD  (NT) 

GB  -  INERTIAL  MAGNETIC  FIELD  GRADIENT  MATRIX  (NT/KM) 

GB(1,1)=DBX/DX 
GB(1,2)=DBX/DY 
GB(1,3)=DBX/DZ 
GB(2,1)=DBY/DX 
GB(2,2)=DBY/DY 
GB(2,3)=DBY/DZ 
GB(3,1)=DBZ/DX 
GB(3,2)=DBZ/DY 
GB(3,3)=DBZ/DZ 

OMEGA  -  MAGNETIC  POTENTIAL  TRANSFORMATION  MATRIX  (KM) 

LAMDA  -  MAGNETIC  FIELD  TRANSFORMATION  MATRIX  (DIMENSIONLESS) 
DLAMD  -  GRADIENT  OF  LAMDA  MATRIX  (KM**(-1)) 


************************************************************************ 


NOTE;  THE  PRISM  IS  ROTATED  THROUGH  EULER  ANGLES  ALPHA, 

BETA  AND  GAMA.  THESE  ANGLES  DEFINE  A  NET 
ROTATION  R(GAMA, BETA, ALPHA) .  THESE  ROTATION  ANGLES 
ARE  DEFINED  IN  ACCORDANCE  WITH  THE  3-2-1  CONVENTION 
THAT  IS  IN  GENERAL  USE  BY  BRITISH  AND  AMERICAN 
AERODYNAMICISTS.  IN  THIS  CONVENTION  THE  ANGLE  ALPHA 
CORRESPONDS  TO  A  COUNTERCLOCKWISE  ROTATION  ABOUT  THE 
POSITIVE  Z-AXIS,  THE  ANGLE  BETA  CORRESPONDS  TO  A 
CLOCKWISE  ROTATION  ABOUT  THE  NEW  Y-AXIS,  AND  THE  ANGLE 
GAMMA  CORRESPONDS  TO  A  COUNTERCLOCKWISE  ROTATION  ABOUT  THE 
FINAL  X  AXIS.  THE  CONSECUTIVE  ROTATIONS  MUST  BE  PERFORMED 
IN  THE  ABOVE  ORDER.  THE  INVERSE  ROTATION  MUST  BE 
PERFORMED  IN  THE  REVERSE  ORDER. 

THE  INERTIAL  COORDINATE  SYSTEM  IS  REFERENCED  TO  AN 
ORIGIN  AT  THE  LOWER  LEFT-HAND  CORNER  OF  THE  SURVEY  AREA. 

THE  PRISM-FIXED  COORDINATE  SYSTEM  IS  REFERENCED  TO  AN 
ORIGIN  THAT  IS  AT  THE  CENTER  OF  THE  PRISM.  THE  PRISM 
FIXED  COORDINATES  ROTATE  WITH  THE  PRISM  RELATIVE  TO 
THE  INERTIAL  COORDINATES. 


********************************************************************* 


SUBROUTINE  MAGREP ( X , Y , Z , V , BX , BY , BZ , GB ) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  MX,MY,MZ,LAMX,LAMY,LAMZ,LAMDA(3, 3) 

DIMENSION  OMEGA (3, 3) , DLAMD (3, 3, 3) ,GB(3,3) , GBP (3, 3) 

DIMENSION  ROT(3,3) ,ROTI(3,3) 

COMMON  /MAGBLK/  MX, MY, MZ, LAMX, LAMY, LAMZ, XB, YB, ZB, ALPHA, BETA, GAMA 
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P1=3.141593D0 
ALPH=ALPHA*PI/ 180 . DO 
BET=BETA*PI/ 180 . DO 
GAM=GAMA*PI/ 180 . DO 

COMPUTE  FORWARD  AND  INVERSE  ROTATION  MATRICES 

ROT(l, l)=+DCOS(ALPH) *DCOS{BET) 

ROT (1,2) =+DSIN ( ALPH) *DCOS (BET) 

ROT(l,3)=-DSIN(BET) 

ROT (2 , 1) =-DSIN (ALPH) *DCOS (GAM) +DCOS (ALPH) *DSIN (BET) *DSIN (GAM) 
ROT (2,2) =+DCOS (ALPH) *DCOS (GAM) +DSIN ( ALPH) *DSIN (BET) *DSIN (GAM) 
ROT (2,3) =+DCOS (BET) *DSIN (GAM) 

ROT ( 3 , 1) =+DSIN (ALPH) *DSIN (GAM) +DCOS (ALPH) *DSIN (BET) *DCOS (GAM) 
ROT (3 , 2 ) =-DCOS (ALPH) *DSIN (GAM) +DSIN (ALPH) *DSIN (BET) *DCOS (GAM) 
ROT(3 , 3 ) =+DCOS (BET) *DCOS (GAM) 

ROTI  ( 1 , 1 )  =+DCOS  ( ALPH )  *  DCOt.  ( BET ) 

ROTI (1,2) =+DCOS (ALPH) *DSIN (BET) *DSIN (GAM) -DSIN (ALPH) *DCOS (GAM) 
ROTI (1, 3) *+DCOS (ALPH) *DSIN(BET) *DCOS (GAM) +DSIN (ALPH) *DSIN(GAM) 
ROTI (2 , 1) =+DSIN (ALPH) *DCOS(BET) 

ROTI (2 , 2 ) =+DSIN (ALPH) *DSIN (BET) *DSIN (GAM) +DCOS (ALPH) *DCOS (GAM) 
ROTI  (2 , 3)  =+DSIN  (ALPH)  *DSIN  (BET)  *DCOS  (GAM)  -DCOS  ( ’'LPH)  *DSIN  (GAM) 
ROTI (3 , 1) *-DSIN (BET) 

ROTI(3,2)*+DCOS(BET) *DSIN(GAM) 

ROTI(3,3)«+DCOS(BET) *DCOS(GAM) 

ROTATE  FROM  INERTIAL  TO  PRISM-FIXED  COORDINATES 

XP»ROT(l,l)*(X-XB)+ROT(l,2)*(y-yB)+ROT(l,3)*(Z-2B) 
yp»ROT (2 , 1) * (X-XB) +ROT (2,2)* (y-yB) +ROT ( 2 , 3 ) * ( Z-ZB) 

ZP*ROT (3,1)* (X-XB) +ROT ( 3 , 2 ) * ( y-yB) +ROT ( 3 , 3 ) * ( Z-ZB ) 

COMPUTE  THE  GEOMAGNETIC  POTENTIAL,  THE  VECTOR  FIELD  AND  THE 
GRADIENT  TENSOR  FIELD  IN  THE  PRISM  FIXED  COORDIANTES 

VP=0.0D0 

BXP=0 . ODO 

ByP=0 . ODO 

BZP^O.ODO 

GBP(1,1)=0.0D0 

GBP(1,2)=0.0D0 

GBP(1,3)=0.0D0 

GBP(2,1)=0.0D0 

GBP(2,2)=0.0D0 

GBP(2,3)=0.0D0 

GBP(3, 1)=0.0D0 

GBP(3,2)=0.0D0 

GBP(3,3)=0.0D0 

DO  30  JB=1,2 

IF  (JB  .EQ.  1)  XL=+LAMX/2.0D0 

IF  (JB  .EQ.  2)  XL=-LAMX/2.0D0 

IF  (JB  .EQ.  1)  SX=+1.0D0 

IF  (JB  .EQ.  2)  SX=-1.0D0 

XD=XP-XL 
DO  20  KB=1,2 

IF  (KB  ,EQ.  1)  yL=+LAMy/2.0D0 

IF  (KB  .EQ.  2)  yL=-LAMy/2.0D0 

IF  (KB  .EQ.  1)  Sy=+1.0D0 

IF  (KB  .EQ.  2)  Sy=-1.0D0 
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YD-YP-YL 
DO  10  LB»1,2 

IP  (LB  .EQ.  1)  ZL=+LAMZ/2.0D0 

IP  (LB  .EQ.  2)  ZL=-LAMZ/2.0D0 

IP  (LB  .EQ.  1)  SZ»+1.0D0 

IP  (LB  .EQ.  2)  SZ=-1.0D0 

ZD«ZP>ZL 
S-SX*SY*SZ 

R2=XD**2+YD**2+ZD**2 

R»DSQRT(R2) 

IP  (R+ZD  .EQ.  O.ODO)  PRINT  *,  ^  A' 

IP  ( (R-YD) / (R+YD)  .EQ.  O.ODO)  PRINT  *,  '  B' 
IP  ( (R-XD) / (R+XD)  .EQ.  O.ODO)  PRINT  *,  '  C' 
IP  (R+XD  .EQ.  O.ODO)  PRINT  '  D' 

LAMDA (1,1) «+S*DATAN ( YD*  ZD/ ( XD*R) ) 
LAMDA(1,2)=-S*DL0G(ABS(R+ZD) ) 

LAMDA (1,3) =-S*DLOG (ABS (R+YD) ) 

LAMDA(2 , 1)  *=LAMDA(  1 , 2) 
LAMDA(2,2)=+S*DATAN(XD*ZD/ (YD*R) ) 

LAMDA (2,3) «-S*DLOG (ABS (R+XD) ) 

LAMDA ( 3 , 1 ) =LAMDA ( 1 , 3 ) 

LAMDA ( 3 , 2 ) *LAMDA ( 2 , 3 ) 

LAMDA (3,3) *+S*DATAN (XD* YD/ ( ZD*R) ) 

OMEGA (1,1) =XD*LAMDA (1,1) 

OMEGA ( 1 , 2 ) «YD*LAMDA ( 1 , 2 ) 

OMEGA (1,3) »ZD*LAMDA (1,3) 

OMEGA (2,1) =XD*LAMDA (2,1) 

OMEGA (2,2) =YD*LAMDA (2,2) 

OMEGA ( 2 , 3 ) =ZD*LAMDA ( 2 , 3 ) 

OMEGA (3,1) «XD*LAMDA (3,1) 

OMEGA ( 3 , 2 ) «YD*LAMDA ( 3 , 2 ) 

OMEGA (3,3) aZD*LAMDA (3,3) 

DLAMD(1,1,2)=-(S/R)*(XD/(R+ZD) ) 

DLAMD(1, 1, 3)=“(S/R) * (XD/ (R+YD) ) 

DLAMD(1,2,1)=DLAMD(1,1,2) 

DLAMD(1,2,2)=-(S/R) *(YD/ (R+ZD) ) 

DLAMD(1,2,3)=-(S/R) 

DLAMD(1,3,1)=DLAMD(1,1,3) 

DLAMD(1,3,2)=DLAMD(1,2,3) 

DLAMD(1,3,3)=-(S/R)*(ZD/(R+YD) ) 

DLAMD (2 , 1, 1) =DLAMD ( 1 , 1 , 2 ) 
DLAMD(2,1,2)=DLAMD(1,2,2) 
DLAMD(2,1,3)=DLAMD(1,2,3) 
DLAMD(2,2,1)=DLAMD(1,2,2) 

DLAMD (2 , 2,3) =- (S/R) * (YD/ (R+XD) ) 
DLAMD(2,3,1)=DLAMD(1,2,3) 
DLAMD(2,3,2)=DLAMD(2,2,3) 
DLAMD(2,3,3)=-(S/R)*(ZD/(R+XD) ) 

DLAMD ( 3 , 1 , 1 ) =DLAMD ( 1 , 1 , 3 ) 
DLAMD(3,1,2)=DLAMD(1,2,3) 
DLAMD(3,1,3)=DLAMD(1,3,3) 
DLAMD(3,2,1)=DLAMD(1,2,3) 

DLAMD(3 , 2 , 2)=DLAMD(2 , 2 , 3) 

DLAMD (3 , 2 , 3) =DLAMD (2 , 3 , 3 ) 
DLAMD(3,3,1)=DLAMD(1,3,3) 
DLAMD(3,3,2)=DLAMD(2, 3,3) 

DLAMD (1,1,1) =-DLAMD (1,3,3) -DLAMD (1,2,2) 
DLAMD(2,2,2)=-DLAMD(2,3,3)-DLAMD(1, 1,2) 
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DLAMD (3,3,3) =-DLAMD (2,2,3) -DLAMD (1,1,3) 

C 

VP=VP-((OMEGA(1,1)+OMEGA(1,2)+OMEGA(1,3) ) *MX 

*  + ( OMEGA (2,1) +OMEGA (2,2) +OMEGA (2,3)) *MY 

*  + ( OMEGA (3,1) +OMEGA (3,2) +OMEGA ( 3 , 3 ) ) *MZ ) 

C 

BXP=BXP+LAMDA (1,1) *MX+LAMDA (1,2) *MY+LAMDA (1,3) *MZ 
BYP*BYP+LAMDA (2,1) *MX+LAMDA (2,2) *MY+LAMDA ( 2 , 3 ) *MZ 
BZP=BZP+LAMDA (3,1) *MX+LAMDA (3,2) *MY+LAMDA (3,3) *MZ 
C 

GBP (1,1) =GBP (1,1) +DLAMD (1,1,1) *MX+DLAMD (1,2,1) *MY+DLAMD (1,3,1) *MZ 
GBP(1,2)=GBP(1,2)+DLAMD(1,1,2)*MX+DLAMD(1,2,2)*MY+DLAMD(1,3,2)*MZ 
GBP(1,3)=GBP(1,3)+DLAMD(1,1,3)*MX+DLAMD(1,2,3)*MY+DLAMD(1,3,3)*MZ 
GBP(2,2)=GBP(2,2)+DLAMD(2,1,2)*MX+DLAMD(2,2,2)*MY+DLAMD(2,3,2)*MZ 
GBP(2,3)=GBP(2,3)+DLAMD(2,1,3)*MX+DLAMD(2,2,3)*MY+DLAMD(2,3,3)*MZ 
GBP(3,3)=GBP(3,3)+DLAMD(3,1,3)*MX+DLAMD(3,2,3)*MY+DLAMD(3,3,3)*MZ 

C 

10  CONTINUE 
20  CONTINUE 
30  CONTINUE 

GBP(2,1)=GBP(1,2) 

GBP(3,1)=GBP(1,3) 

GBP(3,2)=GBP(2,3) 

ROTATE  MAGNETIC  FIELD  COMPONENTS  FROM  PRISM  FIXED  TO 
INERTIAL  COORDINATES 

V=VP 


BX=ROTI (1,1) *BXP+ROTI (1,2) *BYP+ROTI ( 1 , 3 ) *BZP 
BY=ROTI (2 , 1) *BXP+ROTI (2 , 2) *BYP+ROTI ( 2 , 3 ) *BZP 
BZ=ROTI (3,1) *BXP+ROTI (3,2) *BYP+ROTI (3,3) *B2P 

GB(l,l)=ROTI(l,l)*GBP(l,l)+ROTI(l,2)*GBP(2,l)+ROTI(l,3)*GBP(3,l) 
GB(l,2)==ROTI(l,l)*GBP(l,2)+ROTI(l,2)*GBP(2,2)+ROTI(l,3)*GBP(3,2) 
GB(l,3)=ROTI(l,l)*GBP(l,3)+ROTI(l,2)*GBP(2,3)+ROTI(l,3)*GBP(3,3) 
GB(2,1)=!R0TI(2,1)*GBP(1,1)+R0TI(2,2)*GBP(2,1)+R0TI(2,3)*GBP(3,1) 
GB(2,2)*R0TI(2,1)*GBP(1,2)+R0TI(2,2)*GBP(2,2)+R0TI(2,3)*GBP(3,2) 
GB(2,3)=ROTI(2, 1) *GBP ( 1, 3) +ROTI (2 , 2) *GBP(2 , 3) +ROTI (2 , 3) *GBP (3 , 3) 
GB(3,1)=R0TI(3,1)*GBP(1,1)+R0TI(3,2)*GBP(2,1)+R0TI(3,3)*GBP(3,1) 
GB(3,2)=ROTI(3,l)*GBP(l,2)+ROTI(3,2) *GBP(2,2)+ROTI(3,3)*GBP(3,2) 
GB(3,3)=R0TI(3,1)*GBP(1,3)+R0TI(3,2)*GBP(2,3)+R0TI(3,3)*GBP(3,3) 

RETURN 

END 
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c************************************************** ************* ******** 

c 

c 

c  SUBROUTINE  GRVREP  (GRAVITY  FIELD  DUE  TO  A  RECTANGULAR  PRISM) 

C 

C 

C********************************************************** ************* 

c 

c 

C  PROGRAMMED  BY:  JOHN  M.  QUINN  AND  DONALD  L.  SHI EL  8/8/86 

C  NAVAL  OCEANOGRAPHIC  OFFICE 

C  STENNIS  SPACE  CENTER,  MS  39522-5001 

C 
C 

c 

c 

c  PURPOSE:  THIS  ROUTINE  COMPUTES  THE  GRAVITY  FIELD  COMPONENTS  DUE 

C  TO  A  RECTANGULAR  PRISM  LOCATED  AT  (XB,YB,ZB)  RELATIVE  TO 

C  SOME  ORIGIN  LOCATED  AT  THE  OCEAN  SURFACE  AT  THE  LOWER 

C  LEFT-HAND  CORNER  OF  THE  SURVEY  AREA.  THE  PRISM  HAS 

C  DIMENSIONS  (LAMX,LAMY,LAMZ) .  ITS  ORIENTATION  IS 

C  DESCRIBED  BY  EULER  ANGLES  (ALPHA, BETA, GAMA) 

C  CORRESPONDING  TO  YAW,  PITCH,  AND  ROLL  ACCORDING  TO  THE 

C  3-2-1  CONVENTION,  RELATIVE  TO  THE  USUAL  GEODETIC 

C  COORDINATES  FOR  WHICH  X=NORTH,  Y=EAST,  AND  Z=DOWN.  THE 

C  PRISM  HAS  UNIFORM  DENSITY  RHO. 

C 

C 

C 

C 

C  REFERENCE:  A  UNIFIED  APPROACH  TO  GEOPOTENTIAL  FIELD  MODELING; 

C  BY:  JOHN  M.  QUINN  AND  DONALD  L.  SHIEL;  NAVAL 

C  OCEANOGRAPHIC  OFFICE  TECHNICAL  REPORT  #  309  (1993) 

C 

C 

C************************* ************************************* ********* 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

C 

c 

c 

c 

c 

c 

c 

c 


PARAMETER  DESCRIPTIONS: 

X  -  INERTIAL  X  (NORTH)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

Y  -  INERTIAL  Y  (EAST)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

Z  -  INERTIAL  Z  (DOWN)  COORDINATE  OF  OBSERVATION  POINT  (KM) 

XB  -  INERTIAL  X  (NORTH)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

YB  -  INERTIAL  Y  (EAST)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

ZB  -  INERTIAL  Z  (DOWN)  COORDINATE  OF  CENTER  OF  PRISM  (KM) 

XP  -  PRISM  FIXED  X-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

YP  -  PRISM  FIXED  Y-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

ZP  -  PRISM  FIXED  Z-AXIS  COORD.  OF  OBSERVATION  POINT  (KM) 

LAMX  -  DIMENSION  OF  PRISM  ALONG  X-AXIS  (KM) 

LAMY  -  DIMENSION  OF  PRISM  ALONG  Y-AXIS  (KM) 

LAMZ  -  DIMENSION  OF  PRISM  ALONG  Z-AXIS  (KM) 

ALPHA  -  YAW  ROTATION  ANGLE  ABOUT  Z-XIS  OF  PRISM  (DEG.) 

BETA  -  PITCH  ROTATION  ANGLE  ABOUT  Y-AXIS  OF  PRISM  (DEG.) 

GAMA  -  ROLL  ROTATION  ANGLE  ABOUT  X-AXIS  OF  PRISM  (DEG.) 

G  -  NEWTON'S  GRAVITATIONAL  CONSTANT  (KM**3/GM-SEC**2) 

RHO  -  DENSITY  OF  PRISM  (GM/KM**3) 

ROT  -  ROTATION  MATRIX  FROM  INERTIAL  TO  PRISM  FIXED  COORDS. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ROTI  -  INVERSE  ROT.  MATRIX  —  PRISM-FIXED  TO  INERTIAL  COORDS. 

V  -  GRAVITATY  POTENTIAL  IN  INERTIAL  COORDINATES  (KM**2/SEC**2) 

GX  -  INERTIAL  X-COMPONENT  OF  OBSERVED  MAGNETIC  FIELD  {(KM/SEC**2) 

GY  -  INERTIAL  Y-COMPONENT  OF  OBSERVED  MAGNETIC  FIELD  (KM/SEC**2) 

GZ  -  INERTIAL  Z-COMPONENT  OF  OBSERVED  MAGNETIC  FIELD  (KM/SEC**2) 

GG  -  INERTIAL  GRAVITY  FIELD  GRADIENT  MATRIX  (SEC** (-2)) 

GG(1,1)=DGX/DX 
GG(1,2)=DGX/DY 
GG(1,3)=DGX/DZ 
GG(2,1)=DGY/DX 
GG(2,2)=DGY/DY 
GG(2,3)=DGY/DZ 
GG(3, 1)=DGZ/DX 
GG(3,2)=DGZ/DY 
GG(3,3)=DGZ/DZ 

GAMMA  -  GRAVITY  POTENTIAL  TRANSFORMATION  MATRIX  (KM**2) 

OMEGA  -  GRAVITY  FIELD  TRANSFORMATION  MATRIX  (KM) 

LAMDA  -  GRAVITY  GRADIENT  TRANSFORMATION  MATRIX  (DIMENSIONLESS) 


c 

c 

C  NOTE:  THE  PRISM  IS  ROTATED  THROUGH  EULER  ANGLES  ALPHA, 

C  BETA  AND  GAMA.  THESE  ANGLES  DEFINE  A  NET 

C  ROTATION  R (GAMA, BETA, ALPHA) .  THESE  ROTATION  ANGLES 

C  ARE  DEFINED  IN  ACCORDANCE  WITH  THE  3-2-1  CONVENTION 

C  THAT  IS  IN  GENERAL  USE  BY  BRITISH  AND  AMERICAN 

C  AERODYNAMICISTS.  IN  THIS  CONVENTION  THE  ANGLE  ALPHA 

C  CORRESPONDS  TO  A  COUNTERCLOCKWISE  ROTATION  ABOUT  THE 

C  POSITIVE  Z-AXIS,  THE  ANGLE  BETA  CORRESPONDS  TO  A 

C  CLOCKWISE  ROTATION  ABOUT  THE  NEW  Y-AXIS,  AND  THE  ANGLE 

C  GAMMA  CORRESPONDS  TO  A  COUNTERCIX5CKWISE  ROTATION  ABOUT  THE 

C  FINAL  X  AXIS.  THE  CONSECUTIVE  ROTATIONS  MUST  BE  PERFORMED 

C  IN  THE  ABOVE  ORDER.  THE  INVERSE  ROTATION  MUST  BE 

C  PERFORMED  IN  THE  REVERSE  ORDER. 

C 

C  THE  INERTIAL  COORDINATE  SYSTEM  IS  REFERENCED  TO  AN 

C  ORIGIN  AT  THE  LOWER  LEFT-HAND  CORNER  OP  THE  SURVEY  AREA. 

C 

C  THE  PRISM-FIXED  COORDINATE  SYSTEM  IS  REFERENCED  TO  AN 

C  ORIGIN  THAT  IS  AT  THE  CENTER  OF  THE  PRISM.  THE  PRISM 

C  FIXED  COORDINATES  ROTATE  WITH  THE  PRISM  RELATIVE  TO 

C  THE  INERTIAL  COORDINATES. 

C 

c*********************************************** ********************** 


c 

c 

SUBROUTINE  GRVREP (X,Y,Z,V,GX,GY,GZ,GG) 

C 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  RHO, LAMX, LAMY, LAMZ, LAMDA(3 , 3) 

DIMENSION  OMEGA (3, 3) , GAMMA (3, 3) ,GG(3,3) ,GGP(3,3) 

DIMENSION  ROT (3,3), ROTI (3,3) 

COMMON  /GRVBLK/  RHO, LAMX, LAMY, LAMZ, XB, YB, ZB, ALPHA, BETA, GAMA 
C 
C 

PI=3.141593D0 
ALPH=ALPHA*  PI / 1 8  0 . DO 
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BET=BETA*PI/ 180 . DO 
GAM=GAMA*PI / 180 . DO 
G»6.6720D-23 

COMPUTE  FORWARD  AND  INVERSE  ROTATION  MATRICES 

ROT (1,1) =+DCOS ( ALPH) *DCOS ( BET) 

R0T(1,2)*+DS1N(ALPH) *DCOS(BET) 

ROT(l,3)=-DSIN(BET) 

ROT(2 , 1) =-DSIN(ALPH) *DCOS (GAM) +DCOS (ALPH) *DSIN (BET) *DSIN (GAM) 
ROT (2,2) =+DCOS (ALPH) *DCOS (GAM) +DSIN (ALPH) *DSIN (BET) *DSIN(GAM) 
ROT ( 2 , 3) =+DCOS (BET) *DSIN (GAM) 

ROT ( 3 , 1) =+DSIN (ALPH) *DSIN (GAM) +DCOS (ALPH) *DSIN (BET) *DCOS (GAM) 
ROT(3 , 2 ) =- DCOS (ALPH) *DSIN (GAM) +DSIN (ALPH) *DSIN (BET) *DCOS (GAM) 
ROT(3 , 3 ) =+DCOS (BET) *DCOS (GAM) 

ROTI (1,1) =+DCOS (ALPHA) *DCOS (BET) 

ROTI (1,2) =+DCOS (ALPH) *DSIN (BET) *DSIN (GAM) -DSIN (ALPH) *DCOS (GAM) 
ROTI ( 1 , 3) =+DCOS (ALPH) *DSIN (BET) *DCOS (GAM) +DSIN (ALPH) *DSIN (GAM) 
ROTI (2 , 1) =+DSIN (ALPH) *DCOS (BET) 

ROTI (2,2) =+DSIN ( ALPH) *DSIN (BET) *DSIN (GAM) +DCOS (ALPH) *DCOS (GAM) 
ROTI (2 , 3 ) =+DSIN (ALPH) *DSIN (BET) *DCOS (GAM) -DCOS (ALPH) *DSIN (GAM) 
ROTI (3,1) =-DSIN (BET) 

ROTI (3,2) =+DCOS ( BET) *DSIN ( GAM) 

ROTI (3,3) =+DCOS (BET) *DCOS (GAM) 

ROTATE  FROM  INERTIAL  TO  PRISM-FIXED  COORDINATES 

XP«ROT(1,1)*(X-XB)+ROT(1,2)*(Y-YB)+ROT(1,3)*(2-ZB) 

YP«ROT(2,l)*(X-XB)+ROT(2,2)*(Y-YB)+ROT(2,3)*(2-ZB) 

ZP«ROT(3,1)*(X-XB)+ROT(3,2)*(Y-YB)+ROT(3,3)*(Z-ZB) 

COMPOTE  THE  GRAVITY  POTENTIAL,  THE  VECTOR  FIELD  AND  THE 
GRADIENT  TENSOR  FIELD  IN  THE  PRISM  FIXED  COORDIANTES 

VP-O.ODO 

GXP»0 . ODO 

GYP»0 . 000 

GZP«0 . ODO 

GGP(1,1)*0.0D0 

GGP(1,2)*0.0D0 

GGP(1,3)=0.0D0 

GGP(2,1)=0.0D0 

GGP(2,2)»0.0D0 

GGP(2,3)=0.0D0 

GGP(3,1)=0.0D0 

GGP(3,2)=0.0D0 

GGP(3,3)®O.ODO 

DO  30  JB*1,2 

IF  (JB  .EQ.  1)  XL=+LAMX/2.0D0 

IF  (JB  .EQ.  2)  XL*-LAMX/2,0D0 

IF  (JB  .EQ.  1)  SX=+1.0D0 

IF  (JB  .EQ.  2)  SX=-1.0D0 

XD-XP-XL 
DO  20  KB»1,2 

IF  (KB  .EQ.  1)  YL=+LAMY/2.0D0 

IF  (KB  .EQ.  2)  YL=-LAMY/2.0D0 

IF  (KB  .EQ.  1)  SY=+1.0D0 

IF  (KB  .EQ.  2)  SY=-1.0D0 

YD-YP-YL 
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DO  10  LB»1,2 

IF  (LB  .EQ.  1)  ZL— »-LAMZ/2.0D0 

IF  (LB  .EQ.  2)  ZL— LAMZ/2.0D0 

IF  (LB  .EQ.  1)  SZ— t-l.ODO 

IF  (LB  .EQ.  2)  SZ— l.ODO 

ZO-ZP-ZL 
S«SX*SY*SZ 

R2-XD**2+YD**2+ZD**2 

R-DSQRT(R2) 

IF  (R+ZD  .EQ.  O.ODO)  PRINT  * ,  •  K* 

IF  ( (R-YD) / (R+YD)  .EQ.  O.ODO)  PRINT  *,  '  B' 

IF  ( (R-XD) / (R+XD)  .EQ.  O.ODO)  PRINT  '  C' 

IF  (R+XD  .EQ.  O.ODO)  PRINT  * ,  • 

LAMDA(1, 1)-+S*DATAN(YD*2D/ (XD*R) ) 

LAMDA (1,2) — S*DLOG ( ABS (R+ZD) ) 

LAMDA(1,3)— S*DLOG(ABS(R+YD) ) 

LAND A ( 2 , 1 ) -LAMD A ( 1 , 2 ) 

LAMDA (2,2) -+DATAN (XD*ZD/ ( YD*R) ) 

LAMDA(2 , 3) — S*DLOG (ABS (R+XD) ) 

LAMDA ( 3 , 1 ) «LAMDA ( 1 , 3 ) 

LAMDA ( 3 , 2 ) »LAMDA ( 2 , 3 ) 

LAMDA(3,3)=+S*DATAN(XD*YD/ (2D*R) ) 

C 

OMEGA (1,1) «XD*LAMDA (1,1) 

OMEGA ( 1 , 2 ) -YD+LAMDA ( 1 , 2 ) 

OMEGA (1,3) >ZD*LAMDA (1,3) 

OMEGA (2,1) >XD*LAMDA (2,1) 

OMEGA ( 2 , 2 ) «YD*LAMDA ( 2 , 2 ) 

OMEGA ( 2 , 3 ) -ZO*LAMOA ( 2 , 3 ) 

OMEGA (3,1) »XD*LAMDA (3,1) 

OMEGA ( 3 , 2 ) =YD*LAMDA ( 3 , 2 ) 

OMEGA ( 3 , 3 ) -ZD+LAMDA ( 3 , 3 ) 

C 

GAMMA (1,1) »XD*OMEGA (1,1) 

GAMMA ( 1 , 2 ) -YD+OMEGA ( 1 , 2 ) 

GAMMA( 1 , 3 ) »ZD*OMEGA( 1 , 3 ) 

GAMMA ( 2 , 1 ) -XD+OMEGA ( 2 , 1 ) 

GAMMA ( 2 , 2 ) »YD*OMEGA ( 2 , 2 ) 

GAMMA ( 2 , 3 ) »ZD*OMEGA ( 2 , 3 ) 

GAMMA (3,1) sXD+OMEGA (3,1) 

GAMMA ( 3 , 2 ) «YD*OMEGA ( 3 , 2 ) 

GAMMA (3,3) »ZD*OMEGA (3,3) 

C 

VP-VP-G*RHO*((GAMMA(l,l)+GAMMA(l,2)+GAMMA(l,3) ) 

*  +(GAMMA(2,1)+GAMMA(2,2)+GAMMA(2,3) ) 

*  +(GAMMA(3,1)+GAMMA(3,2)+GAMMA(3,3) ) ) 

C 

GXP>=GXP+G*RHO*  (OMEGA  (1,1)  +OMEGA  (1,2)  +OMEGA  (1,3)) 
GYP=GYP+G*RHO* (OMEGA (2,1) +OMEGA (2,2) +OMEGA (2,3)) 
GZP*GZP+G*RHO* (OMEGA (3,1) +OMEGA (3 , 2 ) +OMEGA (3,3)) 
C 

GGP(l,l)=GGP(l,l)+G*RHO*LAMDA(l,l) 

GGP(l,2)»GGP(l,2)+G*RHO*LAMDA(l,2) 

GGP(l,3)=GGP(l,3)+G*RHO*LAMDA(l,3) 

GGP(2,2)»GGP(2,2)+G*RHO*LAMDA(2,2) 

GGP(2,3)»GGP(2,3)+G*RHO*LAMDA(2,3) 

GGP(3,3)=GGP(3,3)+G*RHO*LAMDA(3,3) 

C 

10  CONTINUE 
20  CONTINUE 
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30  CONTINUE 
C 

GGP(2,1)<-GGP(1,2) 

GGP(3,1)-GGP(1,3) 

GGP(3,2)»GGP(2,3) 

ROTATE  GRAVITY  FIELD  COMPONENTS  FROM  PRISM  FIXED  TO 
INERTIAL  COORDINATES 

V-VP 
C 

GX-ROTI (1,1) *GXP+ROTI (1,2) *GYP+ROTI ( 1 , 3 ) *GZP 
GY-ROTI (2 , 1) *GXP+ROTI (2 , 2) *GYP+ROTI (2 , 3 ) *GZP 
GZ-ROTI (3,1) *GXP+ROTI (3,2) *GYP+ROTI ( 3 , 3 ) *GZP 
C 

GG(l,l)»ROTI(l,l)*GGP(l,l)+ROTI(l,2)*GGP(2,l)+ROTI(l,3)*GGP(3,l) 
GG(l,2)*ROTI(l,l) *GGP(l,2)+ROTI(l,2) *GGP(2,2)+ROTI(l,3) *GGP(3,2) 
GG(1,3)«R0TI(1,1)*GGPU»3)+R0TI(1,2)*GGP(2,3)+R0TI(1,3)*GGP(3,3) 
GG(2,l)-ROTI(2,l)*GGP(l,l)+ROTI(2,2)*GGP(2,l)+ROTI(2,3)*GGP(3,l) 
GG(2,2)»R0TI(2,1)*GGP(1,2)+R0TI(2,2)*GGP(2,2)+R0TI(2,3)*GGP(3,2) 
GG(2,3)»ROTI(2,l)*GGP(l,3)+ROTI(2,2)*GGP(2,3)+ROTI(2,3)*GGP(3,3) 
GG(3,1)-ROTI(3,1)*GGP(1,1)+ROTI(3,2)*GGP(2,1)+ROTI(3,3)*GGP(3,1) 
GG(3,2)»ROTI(3,l)*GGP(l,2)+ROTI(3,2)*GGP(2,2)+ROTI(3,3)*GGP(3,2) 
GG(3,3)-R0TI(3,1)*GGP(1,3)+R0TI(3,2)*GGP(2,3)+R0TI(3,3)*GGP(3,3) 
C 

RETURN 

END 
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